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stress  tensor 
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generic  variable  for  illustration  of  finite  volume  integration 
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A  fixed-grid  direct  numerical  simulation  method  has  been  developed  and  verified 
for  single  bubble  dynamics,  heat  transfer  and  phase  change  phenomena.  The  original 
contribution  of  this  research  is  the  first  introduction  of  a  sharp  interface  between  the 
vapor  and  liquid  phases,  which  allows  the  simulation  of  density  ratios  up  to  the  order 
of  1000.  The  interface  location  is  obtained  as  part  of  the  solution  which  accounts  for 
all  coupled  dynamic  and  thermal  effects. 

The  current  method  is  believed  to  be  the  only  method  available  at  present  which 
resolves  the  curvature  calculation  issue,  treats  the  interface  as  a  sharp  discontinuity, 
honors  the  mass  conservation  principle,  is  capable  of  handling  the  large  property 
ratios  and  phase  change. 

For  the  isothermal  cases,  the  predictions  of  drag  coefficients  and  the  shape  defor- 
mation of  bubbles  for  Reynolds  numbers  ranging  from  2  to  100  and  Weber  numbers 
from  2  to  20  were  shown  independently  to  agree  with  previous  published  results. 
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It  was  also  demonstrated  that  the  current  method  is  capable  of  predicting  accu- 
rate results  for  a  wide  range  of  Re,  We,  Ja,  Pe  numbers  and  property  jumps  at  the 
interface. 

For  a  stationary  and  growing  bubble,  the  predicted  growth  rate  was  found  to 
agree  with  the  theoretical  limit  of  oc  t1/2.  However,  for  a  rising  and  growing  bub- 
ble, the  predicted  growth  rate  was  approaching  the  theoretical  limit  of  oc  £2/3  for  a 
spherical  bubble  but  did  not  reach  it  exactly  owing  to  bubble  deformation. 

The  empirical  correlation  proposed  for  bubble  collapse  is  assessed  by  the  simu- 
lations in  the  present  work.  The  present  numerical  model  predicts  that  the  dynamics 
of  a  collapsing  bubble  is  qualitatively,  but  not  quantitatively  similar  to  the  empirical 
correlation  proposed  by  other  investigators. 


CHAPTER  1 
INTRODUCTION 

1.1     Problem  Statement 
1.1.1     Research  Motivation 

Liquid-vapor  phase  change  is  a  phenomenon  that  is  frequently  encountered  in 
daily  life  and  that  plays  an  important  role  in  the  power,  chemical,  petroleum  and 
electronics  industries  among  others.  For  example,  the  power  generation  industry 
takes  advantage  of  the  high  heat  transfer  rates  associated  with  the  liquid- vapor  phase 
change  in  the  heterogeneous  nucleate  boiling  process  to  efficiently  extract  energy 
stored  in  fossil,  nuclear  or  solar  fuels.  The  central  mechanism  of  heat  transfer  during 
nucleate  boiling  is  the  so-called  ebullition  cycle:  a  complete  process  of  liquid  heating, 
nucleation,  bubble  growth  and  departure  (see  Carey,  1992).  Hence  understanding 
and  prediction  of  the  microscopic  vapor  bubble  behavior  have  direct  importance  for 
better  applications  of  nucleate  boiling  heat  transfer. 

Despite  the  wide  applications  of  liquid-vapor  phase  change,  a  fundamental  un- 
derstanding of  such  phase  change  process  is  far  from  being  complete.  Analytical 
analysis  of  the  liquid- vapor  phase  change  has  faced  mathematical  difficulties  because 
of  the  nature  of  nonlinear,  partial  differential  equations  which  govern  the  process. 

In  scientific  and  engineering  research,  experimentation  and  numerical  simulation 
have  been  mainly  employed  as  basic  methods  of  approach.  Both  means  have  respec- 
tive merits  and  they  need  to  complement  each  other.  Meaningful  experiments  are  of 
singular  significance  for  phenomena  that  are  beyond  the  capability  of  current  compu- 
tational technology,  and  they  are  worthwhile  in  providing  verification  for  numerical 


computation  and  guidance  for  mathematical  modeling.  On  the  other  hand,  numeri- 
cal computations  can  be  used  to  predict  microscopic  transport  mechanisms  that  are 
difficult  to  measure  and  visualize  in  an  experiment.  Because  of  limited  resources, 
most  of  the  experimental  projects  may  be  restricted  to  only  one  working  fluid,  simple 
geometry,  and  limited  thermodynamic  conditions.  Numerical  computations  might 
therefore  be  useful  for  expanding  the  data  base  to  a  wider  spectrum  and  providing 
the  information  that  is  otherwise  hard  to  obtain. 

Although  direct  simulations  of  industrial  scale  phase  change  processes  are  cur- 
rently out  of  reach,  a  fundamental  understanding  of  the  complex  physics  at  small 
scale  can  provide  much  insight  to  quantitative  predictions  of  large  scale  aspects  of 
processes. 

Despite  the  fact  that  computational  fluid  dynamics  has  advanced  to  a  point 
where  solutions  of  relatively  complex  flows  without  phase  change  can  be  routinely 
obtained,  numerical  techniques  for  the  combination  of  fluid  flow,  moving  interface 
and  phase  change  remain  as  one  of  the  unresolved  challenges  since  it  is  necessary  to 
model  the  interaction  of  fluid  flow  with  heat  and  mass  transfer  across  moving  phase 
interfaces. 

In  the  present  research,  our  major  goal  is  thus  to  develop  a  numerical  capability 
for  simulating  liquid-vapor  phase  change  at  the  small  scale,  namely,  the  vapor  bubble 
behavior  in  a  pool  of  liquid. 

1.1.2     Physics  of  a  Boiling  Vapor  Bubble  in  Two-Phase  Flows 

Nucleate  boiling  is  one  of  the  most  efficient  heat  transfer  modes.  This  mode  of 
thermal  energy  transport  is  accomplished  by  two  integral  stages,  which  transport  the 
thermal  energy  from  a  heater  surface  to  the  bulk  of  a  working  fluid.  The  first  stage  is 
the  nucleation  and  growth  of  bubble  embryos  on  the  heater  surface.  The  second  stage 
is  the  migration  of  the  bubbles  into  the  bulk  of  the  fluid.  The  present  research  is 
devoted  to  the  fundamental  understanding  of  the  second  stage. 


In  the  second  stage,  if  the  bulk  of  the  liquid  is  superheated,  the  bubble  will 
continue  to  grow  while  it  is  translating.  The  growth  rate  of  the  bubble  in  this  stage 
is  primarily  limited  by  the  rate  of  arrival  of  the  heat  necessary  for  vaporization  at 
the  bubble  surface.  This  stage  of  growth  is  called  heat  transfer-controlled  growth. 
The  pressure  at  the  bubble  surface  inside  the  bubble  approaches  the  liquid  pressure 
adjacent  to  the  bubble  surface;  hence  the  temperature  at  the  bubble  surface  is  ap- 
proximately equal  to  the  equilibrium  saturation  temperature  of  the  system  pressure. 

Translating  vapor  bubbles  in  a  superheated  or  subcooled  liquid  environment  can 
be  found  in  both  pool  boiling  and  convective  boiling.  In  pool  boiling,  bubbles  rise 
into  the  bulk  of  the  pool  due  to  buoyancy.  In  convective  boiling,  bubbles  are  carried 
into  the  bulk  by  flow  currents.  The  interaction  between  the  vapor  bubble  and  its 
liquid  surroundings  takes  place  after  the  departure  of  this  bubble  from  the  heater 
surface.  The  interaction  between  the  dispersed  phase  (bubbles)  and  the  continuous 
phase  (liquid)  involves  exchange  of  momentum,  thermal  energy  and  mass.  The  in- 
teractions are  two-way  coupled  which  means  that  the  bubble  behavior  is  affected  by 
the  liquid  flow  and  the  liquid  flow  is  also  modulated  by  the  presence  of  the  bubbles. 
These  two-way  coupled  interactions  include  momentum  exchange  through  the  drag 
and  vaporization/condensation,  energy  exchange  due  to  convection  and  latent  heat 
transport  of  phase  change,  and  mass  exchange  during  phase  change. 

The  present  research  will  focus  on  numerically  simulating  the  single  vapor  bubble 
growing  or  collapsing  due  to  phase  change  while  translating  in  a  liquid. 

1.1.3     Schematic  of  Computational  Model 

The  schematic  of  computational  model  for  a  growing  or  collapsing  vapor  bubble 
rising  in  a  quiescent  liquid  is  illustrated  in  Figure  1.1.  The  motion  of  the  vapor 
bubble  is  simulated  in  a  cylindrical  domain  filled  with  liquid  of  the  same  substance. 
The  liquid  is  superheated  in  one  case  or  subcooled  in  the  other  case.  The  major 
simplification  is  that  the  system  is  assumed  axisymmetric. 


Figure  1.1:  Schematic  of  computational  model. 


1.1.4     Assumptions 

The  two-fluids  model  is  adopted  as  the  mathematical  basis  which  governs  the 
liquid- vapor  two-phase  flow  of  a  single  translating  vapor  bubble;  i.e.,  a  separate  set 
of  mass,  momentum  and  energy  conservation  equations  are  solved  for  each  phase 


(dispersed  vapor  phase  and  continuous  liquid  phase  surrounding  the  bubble)  while 
sharp  discontinuities  of  material  properties  are  maintained  across  the  interface  of  zero 
thickness. 

The  assumptions  introduced  throughout  this  study  are  summarized  as  follows: 

•  Axisymmetric,  cylindrical  coordinate  system 

•  Single  component  system 

•  Newtonian,  constant  property  fluids  in  each  phase 

•  Thermal  equilibrium  at  the  phase  interface  which  means  temperatures  at  the 
interface  in  the  liquid  and  vapor  phases  are  equal 

1.1.5     Governing  Equations 

With  the  above  assumptions,  the  mass,  momentum  and  energy  conservation 
equations  for  both  phases  can  be  written  as 

V-u   *   0  (l.l) 

du 


+  V-(««) 


=    -Vp  +  V  ■  (/iVu)  (1.2) 


pcp 


f+v.(«D 


=    V-(fcVT)  (1.3) 


where  for  numerical  efficiency,  pg  =  -A(pgh)  has  been  combined  with  p  in  place 
of  the  pressure.  If  the  actual  pressure  is  needed,  hydrostatic  pressure  pgh  has  to  be 
subtracted  from  the  pressure  computed  in  the  equations.  Thereby  the  hydrostatic 
pressure  (body  force  term)  contribution  will  only  appear  in  the  normal  stress  balance 
at  the  bubble  interface  (see  Eq.  (1.12)). 

For  this  research,  the  axisymmetric  cylindrical  coordinates  are  employed.  Equa- 
tions (1.1)-(1.3)  thus  become  the  following  with  d/d6  =  0  and  the  swirl  velocity 
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where  the  stress  tensor  components  are 
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1.1.6     Boundary  and  Interfacial  Conditions 


Apart  from  the  governing  equations,  boundary  conditions  are  needed  to  make 
the  problem  well  posed.  The  interfacial  conditions,  which  are  simply  the  statement 
of  those  conservation  principles  at  the  interface,  must  also  be  satisfied  (see,  e.g., 
Delhaye,  1974;  Ishii,  1975;  Kataoka,  1986;  Sadhal,  Ayyaswamy  &  Chung,  1996). 
Boundary  conditions: 

The  outflow  (zero  gradient  of  the  velocity)  conditions  are  specified  on  the  top, 
left  and  right  sides  of  the  domain.  At  the  bottom  boundary,  the  symmetric  condition 
for  all  variables  is  used. 


Interfacial  mass  conservation  conditions: 

Due  to  phase  change,  for  an  interface  moving  at  a  velocity  Uint,  the  mass  con- 
servation conditions  on  the  interface  stipulates  that 

rh  =  pi  (uint  -  ui)  n    =    pv  (uint  -  uv)  n 
=>    rh  =  pi  [(un)int  -  (un)i]    =    pv  [(un)int  -  («„)„]  (1.9) 

where  rh  is  the  interfacial  mass  flux  (positive  for  evaporation),  pt  and  pv  are  the  liquid 
and  vapor  densities  respectively,  while  the  tangential  velocity  is  always  continuous 

urt  =  uvt    =s»    (ut)i  =  (ut)v  (1.10) 

Interfacial  momentum  balance  conditions: 

The  forces  balance  across  the  interface  is  expressed  as 

Pv  [{Un)v  -  Mint]  Uv  +  PVn  -  TVU     =     pt  [(un)i  -  (un)int\  Ut 


+PlTl  —  TrU  +  OKU 


(1.11) 


The  normal  component  of  the  equation  (1.11)  is 

Pi  -  pv  +  an    =    (Trn)n  -  (rv-n)n  -  {pt  [(un)i  -  (un)int]  (un)t  - 

Pv  [(Un)v  ~  Mint]  K)„}  (1.12) 

where  r  represents  the  viscous  stress  tensor,  a  is  the  surface  tension  coefficient  which 
is  assumed  to  be  constant,  k  =  (-^  +  j^)  with  i?i  and  R2  the  principal  radii  of 
interface  curvature.  The  pressure  in  Eq.  (1.12)  is  the  total  pressure  including  the 
hydrostatic  contribution;  i.e.,  the  actual  pressure  is  obtained  by  subtracting  pgh  (or 
j^h  in  dimensionless  form)  from  p  solved  in  Eq.  (1.1),  (1.2). 

The  normal  component  of  stress  balance  is  reduced  to  Young-Laplace  equation 
when  there  is  no  motion  in  both  phases. 
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Interfacial  energy  conservation  conditions: 

The  energy  conservation  relation  across  the  interface  involves  both  the  sensible 
and  latent  heat  portions 

-  (kvVTv  -  hVTi)  n    =    -Xpv  (uv  -  uint)  n 

fdT\  (dT\ 

SSS>    kl\~d~)    ~kv\d~)      =    ^Pv[{u>n)int-{un)v\  =  \rn  (1.13) 

where  A  is  the  latent  heat  of  evaporation. 

Equation  (1.13)  is  used  to  estimate  the  normal  velocity  component  of  the  inter- 
face driven  by  the  phase  change. 

Interface  temperature  conditions: 

By  thermal  equilibrium,  the  temperature  continuity  condition  reads  as 

Tt  =  Tv  =  Tint  (1.14) 

where  Tint  is  the  temperature  at  the  phase  change  interface.    A  condition  on  the 
interface  temperature  must  be  specified  to  complete  the  formulation. 

In  the  case  where  liquid  boiling  occurs  on  a  horizontal  flat  surface,  i.e.,  curvature 
k  —  0,  the  interface  temperature  Tint  is  equal  to  saturation  temperature  Tsat  of 
the  two  phase  mixture  at  the  corresponding  ambient  pressure  p^.  For  a  curved 
interface  between  the  liquid  and  vapor  bubbles,  when  surface  tension  effect  is  taken 
into  account,  the  thermal  and  mechanical  equilibrium  at  the  interface  reveals  that  the 
local  saturation  temperature  Tint  along  the  interface  will  deviate  slightly  from  that 
constant  saturation  temperature  Tsat  for  the  flat  interface  (see  Alexiades  &  Solomon, 
1993),  which  is  the  so-called  Gibbs-Thomson  effect.  For  the  heat  transfer  controlled 
bubble  growth,  vapor  pressure  pv  at  the  bubble  interface  is  equal  to  the  system 
pressure  p^  which  is  the  equilibrium  pressure  of  Tsat,  i.e., 

Pv=Poo  (1.15) 


In  this  case,  a  sufficiently  accurate  formulation  for  the  interface  phase  change  tem- 
perature reads 

Tint  =  Tsat(l  +  ^)  (1.16) 

1.1.7     Nondimensionalization 

To  nondimensionalize  the  governing  equations  in  the  liquid  and  vapor  phases 

as  well  as  boundary  and  interfacial  conditions,  the  reference  scales  are  length  L, 

velocity  Ur,  time  tT  —  L/uT.  The  characteristic  temperature  scale  is  AT  —  T)l  —  Tc  = 
Too  ~  Tsat    growth 


( 


.    Too  and  Tsat  are,  respectively,  ambient  temperature  and 
Tsat -Too    collapse 

saturation  temperature. 

With  the  basic  scales  chosen  as  above,  the  dimensionless  variables  are  defined 
as  x*  =  x/L,  u*  =  u/ur,  t*  =  t/tr,  p*  =  p/{piur2),  T*  =  (T  -  TC)/AT, 
9*  =  9/9,  P*  =  P/Pref,  P*  =  p/Pref,  k*  =  k/kref,  cp*  =  cp/cPref,  where  the 
liquid  properties  are  taken  to  be  the  reference  properties,  i.e.,  pTej  =  pi,  p,Tef  =  pi, 

™ref  ~~  "<Jj  Cpref  —  Cpi- 

Finally,  the  dimensionless  equations  in  the  liquid  and  vapor  phases,  omitting  the 
superscript  *,  can  be  written  as 

Liquid  Phase: 

V-u    =    0  (1.17) 

-57  +  V-(ttt»)    =    -Vp+—  V2u  (1.18) 

at  Re 

dT  1 

—  +  V-(«T)    =    -V2T  (1.19) 

Vapor  Phase: 

V-u    =    0  (1.20) 


Py_ 
.Pi. 


du 

lit 


+  V  •  (uu) 


MS*;*" 
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(1.21) 


dT 
^+V.(,T) 


k, 


■Pi 
"-Vv 


c„..  I  Pe 


±V2T 


where 


Reynolds  number        Re 


Peclet  number         Pe 


pturL 

Pi 
picPlurL 


Interfacial  Condition: 

Mass  continuity  condition 

\Un)int      = 

Momentum  balance  condition 


Mi  ~  (g)  M„ 


n-K  +  We* 


1 


dun 
dn  }x 


pA    fdUn 

fn)  \dn 

'  Pv^ 


(1.22) 


(1.23) 


[Ml  -  Mint)  Ml  +  (  —  J  [M„  -  Mint)  M„ 


(1.24) 

where  pi  =  Pd  +  If  Ft  •  g{zout  -  z),    Fr  =  u2TjgL  is  the  total  pressure  including 
hydrostatic  and  dynamics  pressure.  zouf  is  the  level  of  the  liquid  pool. 
Energy  conservation  condition 


\Un)int        \Un)v  ~  JT  ' 

Interface  temperature  condition 


37}       (kv\  (dTv 


dn       \ki  J  \dn 


(1.25) 


J-int  —  '  1$ 


(1.26) 
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where 


Weber  number  We  = 

Froude  number  Fr  = 

Jakob  number  Ja  = 

Parameter  J1  = 


Piu;L 


a 

gL 

pi 

cPlAT 

Pv 

T 

±  sat 

A 
a 

AT  PlXL 


1.1.8     Dimensionless  Parameters 


After  examining  the  dimensionless  governing  equations  and  the  associated  bound- 
ary and  interfacial  conditions,  the  major  dimensionless  parameters,  namely,  Re,  Fr, 
Pe,  Pr,  We  and  Ja,  and  their  respective  physical  meanings  are  listed  as  follows: 

piurL 
Re    =    ,         Reynolds  number 


flow  inertia  force 

viscous  force 
i 

'7,         Froude  n 

J 

flow  inertia  force 


v? 
Fr    —    — £-,         Froude  number 
9L 


gravity  force 

Pe    —    — m  T    ,         Peclet  number 
ki 
convection 


conduction 
P\urL    yncVl 

Pi         h 

RePr 


Pr    =    — ,      vi  =  — ,     ax  = 1         Prandtl  number 

on  Pi  PicPl 

momentum  diffusion  capacity 
thermal  diffusion  capacity 

We    =    — — ,         Weber  number 
a 

flow  hydrodynamic  pressure  force 
surface  tension  force 
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J  a    =    —  •  -3EU — ,         Jakob  number 
Pv        A 

sensible  heat  transfer 
latent  heat  transfer 

1.2     Literature  Review 
1.2.1     Vapor  Bubble  Growth  or  Collapse  in  a  Two-Phase  Flow 

Since  we  are  concerned  with  the  heat  transfer-controlled  growth  or  collapse  of 
a  single  vapor  bubble  translating  in  a  liquid  pool,  the  review  of  previous  works  will 
be  focused  only  on  those  directly  relevant  to  thermally  controlled  bubble  growth  or 
collapse  after  the  bubble  has  departed  from  the  heater  surface.  A  number  of  excellent 
publications  in  the  literature  are  dedicated  entirely  to  the  fundamental  understand- 
ing of  the  bubble  growth  dynamics  on  the  heater  surface  (first  stage  of  boiling),  for 
example,  Klausner,  Mei  &  Bernhard  (1993);  Zeng,  Klausner  &  Mei  (19936);  Zeng, 
Klausner,  Bernhard  &  Mei  (1993a);  Chen,  Klausner  &  Mei  (1995)  and  Mei,  Chen  k 
Klausner  (1995).  These  works  are  not  directly  relevant  to  the  scope  of  the  proposed 
research,  and  will  not  be  discussed  in  detail  further  here. 

Vapor  bubble  growth  or  collapse  has  drawn  extensive  study  over  many  years 
because  of  its  theoretical  implication  to  the  nucleate  boiling  and  other  two  phase 
physical  phenomena.  Due  to  its  intrinsic  difficulties,  mainly  the  nonlinearity  of  con- 
vection and  the  coupling  of  the  bubble  evolution  with  the  surrounding  liquid  motion 
as  well  as  the  large  number  of  variables  involved,  it  would  seem  that  this  problem 
has  been  resolved  to  a  tenable  degree  of  understanding. 

The  heat  transfer  controlled  bubble  growth  or  collapse  can  be  generally  divided 
into  four  basic  situations:  stationary  bubble  growth,  translating  bubble  growth,  sta- 
tionary bubble  collapse  and  translating  bubble  collapse.  The  researchers  in  this  field 
seem  to  follow  this  line  to  report  results  for  one  or  more  cases  in  any  particular  paper. 
In  our  discussion,  we  will  also  follow  this  convention  to  look  at  the  prior  contributions 
to  all  these  cases. 
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Stationary  bubble  growth 

The  majority  of  the  research  on  the  heat  transfer  controlled  growth  of  station- 
ary vapor  bubbles  is  based  on  the  simplified  model  that  a  spherical  bubble  in  an 
unbounded  fluid  undergoes  spherical  growth.  The  simplified  model  does  not  solve  for 
the  hydrodynamics  of  the  fluid  flow,  nor  does  it  locate  the  moving  phase  interface. 
With  the  assumption  of  spherical  symmetry,  the  continuity  and  momentum  equa- 
tions for  an  incompressible  fluid  lead  to  the  famous  Rayleigh  equation  (see  Plesset  & 
Prosperetti,  1977) 

RW+i\&)    =«ift(T«)-p--B-4flUJ}  (127) 

in  which  Ts(t)  is  the  liquid  temperature  at  the  bubble  surface  and  it  is  assumed  that 

the  vapor  pressure  is  the  thermodynamic  equilibrium  pressure  at  the  temperature  Ts. 

To  analyze  the  heat  transfer  controlled  growth,  the  liquid  temperature  Ts  at  the 

bubble  surface  is  obtained  by  solving  the  energy  equation 

dT  ,  R2dRdT      a,  d  (  2dT 


dt       r2  dt  dr       r2  dr  \     dr  J 
in  which  ai  —  ki/(picpi)  is  the  liquid  thermal  diffusivity. 

Plesset  &  Zwick  (1954),  Forster  &  Zuber  (1954),  Birkhoff,  Margulies  &  Horning 
(1958)  and  Scriven  (1959)  considered  the  heat  transferred  controlled  spherical  vapor 
bubble  growth  without  translation  by  solving  the  energy  equation  (1.28)  subject 
to  certain  approximations  in  conjunction  with  the  dynamical  equation  (1.27)  and 
thermodynamic  equilibrium  relation  pv  =  pv(Ts).  Their  results  of  bubble  radius 
varying  with  time  are  in  good  agreement  with  the  experimental  data  of  Dergarabedian 
(1953)  and,  later,  Florschuetz,  Henry  &  Khan  (1969);  basically,  what  they  concluded 
is  that  the  bubble  radius  of  the  heat  transfer  controlled  growth  is  proportional  to 
(time)1'2. 

Prosperetti  &  Plesset  (1978)  re-examined  two  assumptions  introduced  in  those 
earlier  analyses  of  bubble  growth:  the  validity  of  a  thin  thermal  boundary  layer  near 
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Figure  1.2:  Bubble  radius  vs.  time  for  a  vapor  bubble  growing  in  superheated  water 
at  atmospheric  pressure  with  5°C  superheat  [Prosperetti  &  Plesset  (1978)]. 


the  bubble  wall  and  the  empirical  linear  relation  for  the  vapor  pressure  with  vapor 
temperature.  It  was  showed  that,  upon  comparison  with  numerical  results  of  Donne 
&  Ferranti  (1975),  both  assumptions  do  not  lead  to  serious  errors  in  the  radius-time 
behavior  for  superheat  ranges  of  practical  interest.  The  relation  of  the  bubble  radius 
with  time  for  heat  transfer  controlled  growth  is  given  as 


R  =  2kt 


sat 


t* 


(1.29) 


_  -non  j     \      Xpv 

in  which  fcj  is  the  liquid  thermal  conductivity,  at  =  ki/(picPl)  is  the  liquid  thermal 
diffusivity,  (Too  -  Tsat)  is  the  liquid  superheat,  A  is  the  latent  heat,  pv  is  the  vapor 
density.  The  radius  increasing  with  time  based  on  this  formulation  is  plotted  in 
Figure  1.2. 
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Translating  bubble  growth 

The  stationary  vapor  bubble  is  an  ideal  exception.  In  practice,  the  vapor  bubble 
usually  exhibits  more  or  less  translatory  motion.  Although  the  translatory  motion  is 
not  expected  to  play  a  significant  role  if  liquid  inertia  is  the  dominant  mechanism, 
as  stated  by  several  investigators,  it  has  a  large  influence  on  the  growth  or  collapse 
if  heat  transfer  is  controlling  the  phase  change.  Under  such  circumstances,  the  heat 
transfer  would  be  enhanced  by  translation,  thus  resulting  in  higher  growth  or  collapse 
rates. 

The  experimental  investigation  of  Darby  (1964)  studied  simultaneous  growth  and 
translation  characteristics  of  vapor  bubbles  issuing  continuously  from  a  nucleation 
site.  The  bubble  radius  was  found  to  vary  as  (time)2/3  during  the  entire  recorded 
history.  It  indicates  that,  as  expected,  the  translational  motion  increases  the  growth 
rate  over  that  for  stationary  bubbles. 

Ruckenstein  &  Davis  (1971)  assumed  a  potential  flow  in  the  region  surrounding 
the  bubble  to  analyze  the  effect  of  bubble  translation  on  vapor  bubble  growth  in  a 
superheated  liquid.  He  found  that  for  a  moderate  Jakob  number,  the  translation  can 
substantially  increase  bubble  growth  rate  over  the  rates  predicted  for  a  stationary 
growing  bubble.  For  sufficiently  large  Jakob  numbers  (J a  >  50),  bubble  growth  is  not 
greatly  affected  by  the  relative  motion  between  the  bubble  center  and  the  surrounding 
liquid.  The  results  predicted  by  Ruckenstein  &  Davis  (1971)  are  in  good  agreement 
with  the  experimental  data  of  Florschuetz,  Henry  k,  Khan  (1969). 

In  Table  1.1  we  summarize  the  above  mentioned  important  research  publications 
for  heat  transfer  controlled  growth  of  a  vapor  bubble. 

Stationary  bubble  collapse 

In  bubble  collapse,  the  radius  variation  produced  by  a  given  interfacial  heat  flux 
increases  in  time  in  contrast  to  the  decrease  for  the  bubble  growth  case.  This  interface 
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Table  1.1:  A  list  of  representative  investigations  on  stationary  and  translating  vapor 
bubble  growth 


Authors 

Approach 

Assumptions  &  Conditions 

Major  Results 

Plesset  &  Zwick 

Analytical 

-  Spherical  symmetry 

The           bubble 

(1954) 

-  Stationary  growing  vapor  bubble 

radius     is     pro- 

Forster &  Zuber 

-  Heat  transfer  controlled  vapor  bub- 

portional         to 

(1954) 

ble  growth 

the   square  root 

Birkhoff,      Mar- 

-  The  inertia  and  viscosity  of  the  liq- 

of     time,       i.e. 

gulies  &  Horning 

uid  are  neglected 

(time)1/2 

(1958) 

-  Surface  tension  is  neglected 

Scriven  (1959) 

-  Heat  transfer  by  conduction  only 

Prosperetti       & 

-    Constant    thermal    and    physical 

Plesset  (1978) 

properties 

Darby  (1964) 

Experiment 

it  Translating  vapor  bubble  growth 

Bubble  ra- 
dius varies  as 
(time)2/3 

Ruckenstein      &: 

Analytical 

-  Spherical  symmetry 

The           growth 

Davis  (1971) 

-  Translating  vapor  bubble  growth 

rate  is  substan- 

- Potential  flow  in  the  liquid 

tially    increased 

-  Heat  transfer  controlled  growth 

at  moderate 
Jakob  number 
{Ja  <  50) 

Legendre,  Boree 

Numerical 

-  Spherical  symmetry 

Confirmed       all 

&       Magnaudet 

-   Full  Navier-Stokes  equations  are 

above  prior  an- 

(1998) 

solved 

alytical     results 

-  Heat  transfer  controlled  growth 

and  studied  the 

-  Stationary  and  translating  bubble 

hydrodynamic 

-   Decoupled   bubble  evolution   and 

forces 

liquid  motion 

behavior  makes  condensation  basically  more  difficult  to  describe  than  the  evapora- 
tion. For  the  stationary  bubble  collapse,  by  virtue  of  spherical  symmetry  assumption, 
Florschuetz  &  Chao  (1965)  examined  the  relative  importance  of  the  liquid  inertia  and 
heat  transfer  effects  on  the  vapor  bubble  collapse.  A  dimensionless  parameter  is  iden- 
tified to  characterize  the  collapse  mode:  either  liquid  inertia-controlled  collapse,  heat 
transfer  dominated  collapse  or  an  intermediate  case  where  both  liquid  inertia  and 
heat  transfer  are  of  comparable  importance.    For  heat  transfer-controlled  collapse, 
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they  found  that  the  collapse  rate  is  relatively  slow  and  decreases  as  the  collapse  pro- 
ceeds, the  radius-time  relation  is  bounded  by  two  limiting  heat  transfer-controlled 
collapse  curves  that  resulted  from  two  approximations  for  temperature  solutions: 

7  =  l-r^  (1.30) 

and 

x„  =  i(?+72-3)  (1.31) 

in  which  dimensionless  bubble  radius  7  =  R/Ro,  dimensionless  time 

rn  =  ija'£  (1.32) 

with  Ja  being  the  Jakob  number,  at  the  liquid  thermal  diffusivity,  Rq  the  initial 
radius  of  the  bubble.  The  results  are  visualized  in  Figure  1.3  which  is  from  the  paper 
of  Florschuetz  &  Chao  (1965). 

The  analytical  results  Eq.  (1.31)  of  Florschuetz  &  Chao  (1965)  are  valid  at  high 
Jakob  numbers.  In  order  to  obtain  results  valid  at  finite  Jakob  numbers,  Okhotsimskii 
(1988)  solved  numerically  the  heat  transfer  controlled  collapse  based  on  the  same 
theoretical  model  of  Plesset  &  Zwick  (1954)  and  tabulated  the  time  evolution  of  the 
bubble  radius  over  a  wide  range  of  Jakob  numbers  (  10~2  <  J  a  <  103  ).  In  his 
analysis,  he  isolated  the  Jakob  number  to  be  the  only  dominant  factor  of  the  collapse 
process. 

Translating  bubble  collapse 

For  more  practical  cases  of  translating  vapor  bubble  collapse,  Wittke  k  Chao 
(1967)  obtained  numerical  and  experimental  information  of  the  influence  of  transla- 
tory  velocity  on  the  collapse  of  a  vapor  bubble.  Noting  that  the  translatory  motion  is 
not  expected  to  be  a  significant  factor  if  liquid  inertia  is  the  dominant  mechanism  for 
collapse,  they  examined  the  bubble  traversing  at  a  constant  velocity  collapsing  under 
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Figure  1.3:  Bubble  radius  vs.  time  for  heat  transfer  controlled  vapor  bubble  collapse 
[Florschuetz  &  Chao  (1965)]. 


heat  transfer  controlled  conditions.  They  found  that,  again  as  expected,  because  the 
relative  contribution  to  heat  transfer  at  the  bubble  surface  due  to  translatory  mo- 
tion is  greater  when  the  Jakob  number  is  smaller,  the  translatory  motion  produces  a 
greater  effect  for  the  smaller  Jakob  number.  Overall,  the  faster  the  bubble  translates, 
the  more  rapidly  it  collapses. 

Tokada,  Yang  &  Clark  (1968)  and  Dimic  (1977)  also  studied  theoretically  the  ef- 
fect of  translational  velocity  of  a  bubble  on  the  collapse  rate.  They  basically  obtained 
similar  results  on  the  collapse  rates. 

In  these  studies,  the  bubble  is  relatively  large  (2  mm  <  Ro  <  4  mm)  and  because 
of  reasons  which  are  not  completely  clear,  it  was  observed  in  some  experiments  on 
bubble  collapse  that  bubbles  exhibit  almost  a  constant  rising  velocity  (see  Wittke 
&  Chao,  1967;  Chen  &  Mayinger,  1992);  hence  a  constant  translation  velocity  is 
assumed  in  all  these  studies. 
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Moalem  &  Sideman  (1973)  analyzed  the  effect  of  variable  rising  velocities  as- 
sociated with  relatively  small  bubbles  (R0  <  1  mm)  as  well  as  the  effect  of  a  cross 
flow;  their  results  were  in  good  agreement  with  the  experimental  data  of  Abdelmes- 
sih,  Hooper  &  Nanngia  (1972).  What  they  concluded  is  that  the  effect  of  variable 
velocity  as  compared  to  a  constant  velocity  motion  is  relatively  small. 

All  of  the  above  mentioned  results  assumed  the  potential  flow  for  outer  liquid 
surrounding  the  bubble. 

Chen  &  Mayinger  (1992)  performed  very  sophisticated  measurements  to  deter- 
mine the  heat  transfer  at  the  phase  interface  of  a  collapsing  vapor  bubble  rising  in 
a  subcooled  liquid.  The  temporal  decrease  in  the  bubble  diameter  is  given  in  the 
correlation  using  experimental  data  as 

7  =  (1  -  0.56Re°-7Pr0-5JaFo)°-9  (1.33) 

in  which  7  =  R/Ro,  Reynolds  number  is  based  on  initial  bubble  diameter.  The 
Fourier  number  is 

F°  =  W^  ^ 

Recently,  Gumerov  (1996)  numerically  obtained  the  solutions  of  the  same  theoretical 
model  as  that  of  Wittke  &  Chao  (1967)  for  the  collapse  of  a  vapor  bubble  with 
constant  translatory  velocity.  His  results  are  shown  to  agree  substantially  better 
with  the  experiments  in  a  number  of  cases. 

In  Table  1.2  some  of  these  studies  for  stationary  and  translating  vapor  bubble 
collapse  are  summarized. 

Most  recently,  Legendre,  Boree  &  Magnaudet  (1998)  made  a  comprehensive 
study  of  heat  transfer-controlled  spherical  bubble  growth  and  collapse  without  and 
with  constant  translation  velocity.  Basically,  they  reexamined  all  four  cases,  and  by 
comparing  with  those  representative  prior  results,  confirmed  previous  analytical  anal- 
ysis. In  their  study,  a  numerical  approach  is  employed  to  solve  the  full  Navier-Stokes 
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Table  1.2:  A  list  of  representative  investigations  on  stationary  and  translating  vapor 
bubble  collapse 


Authors 

Approach 

Assumptions  &  Conditions 

Major  Results 

Florschuetz  & 

Analytical 

-  Spherical  symmetry 

Three  collapse 

Chao  (1965) 

and 

-  Stationary  collapsing  bubble 

modes          are 

Experimental 

-  Surface  tension  and  viscosity  are  ig- 

classified and  a 

nored 

dimensionless 

-  No  relative  velocity  between  liquid 

parameter       is 

and  bubble  wall 

identified       to 

-  Heat  conduction  in  vapor  is  ignored 

characterize 

-  Pressure  in  the  vapor  phase  is  uni- 

the       collapse 

form 

mode. 

-  Perfect  gas  for  vapor 

Okhotsimskii 

Numerical 

-  Spherical  symmetry 

Time     evolu- 

(1988) 

-  Stationary  collapsing  vapor  bubble 

tion  of  bubble 

-  Heat  transfer  controlled  vapor  bub- 

radius   for    a 

ble  collapse 

wide        range 

-  The  inertia  and  viscosity  of  the  liq- 

of          Jakob 

uid  are  neglected 

numbers        is 

-  Surface  tension  is  neglected 

tabulated. 

-    Constant    thermal    and    physical 

properties 

Wittke          & 

Numerical 

-  Spherical  symmetry 

The       transla- 

Chao (1967) 

and 

-  Translating  vapor  bubble 

tion        motion 

Experimental 

-  Heat  transfer  controlled  collapse 

generally       in- 

- Constant  translational  velocity 

creases         the 

-  Potential  outer  flow 

collapse  rates. 

-  Incompressible  for  Liquid 

Chen             & 

Experimental 

-  Translating  vapor  bubble 

Correlations  of 

Mayinger 

-  Heat  transfer  controlled  collapse 

Nusselt     num- 

(1992) 

ber    and    time 
evolution       of 
bubble    radius 
are  given. 

Legendre, 

Numerical 

-  Spherical  symmetry 

Confirmed     all 

Boree            & 

-  Full  Navier-Stokes  equations  are 

above        prior 

Magnaudet 

solved 

analytical 

(1998) 

-  Heat  transfer  controlled  collapse 

results         and 

-  Stationary  and  translating  bubble 

studied        the 

-  Decoupled  bubble  evolution  and 

hydrodynamic 

liquid  motion 

forces 
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Figure  1.4:  Computational  configuration  of  investigation  of  Legendre,  Boree  &  Mag- 
naudet  (1998). 


equations  and  temperature  equation  in  the  liquid  to  determine  the  temperature  field 
around  the  bubble  surface  so  as  to  obtain  the  growth  or  collapse  rates.  Despite  the 
spherical  symmetry  assumption  and  disregarding  of  the  coupling  between  the  bubble 
evolution  and  the  liquid  motion,  their  solutions,  to  our  knowledge,  are  the  most  ad- 
vanced ones  for  the  spherical  bubble  growth  and  collapse.  The  second  contribution 
of  their  research  is  regarding  the  hydrodynamic  force  experienced  by  the  bubble  with 
a  time-dependent  radius. 

The  computational  configuration  of  Legendre,  Boree  &:  Magnaudet  (1998)  is 
shown  in  Figure  1.4. 
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1.2.2     Numerical  Methods  for  Moving  Interface  Problems 

The  liquid-vapor  two-phase  problem  we  outlined  in  chapter  1  is  a  special  case 
of  more  general  situations  involving  multi-fluids  or  multi-phases  in  a  computational 
domain  as  shown  in  Figure  1.5,  where  multiple  immersed  interfaces  are  embedded 
in  the  flow  field.  The  regions  enclosed  by  the  interfaces  designate  certain  phases 
which  could  be  solid,  liquid  or  gas  phases  and  might  be  different  from  that  of  the 
surrounding  fluid.  The  interfaces  can  be  stationary,  or  moving,  deforming  through 
interaction  with  the  flow  around  it.  Such  a  scenario  may  be  encountered  in  several 
types  of  flow  problems.  For  instance,  the  interface  may  represent  a  solid  geometry 
over  which  fluid  flows.  This  solid  boundary  may  be  rigid  as  in  pure  external  flow 
problems  or  be  flexible  exhibiting  motions  under  the  influence  of  the  external  flow  as 
in  many  fluid-structure  interaction  problems.  In  other  situations,  the  interface  may 
represent  a  phase  boundary  at  which  a  phase  transition  occurs.  One  of  the  examples 
is  the  liquid-vapor  phase  change  problem  involving  vapor  bubbles  as  in  the  present 
research,  the  enclosed  phase  is  the  dispersed  vapor  phase,  and  the  surrounding  one 
is  the  continuous  liquid  phase.  In  an  inverse  situation  involving  liquid  drops,  the 
enclosed  phase  is  the  liquid  phase,  and  continuous  phase  is  the  vapor  phase.  In 
these  cases,  sharp  discontinuities  of  the  fluid  properties  such  as  density,  viscosity  and 
thermal  conductivity  arise  across  the  interface.  Other  examples  of  phase  boundary 
include  solid-liquid  phase  of  the  solidification  in  material  processing  problems.  The 
interface  may  also  simply  represent  a  boundary  separating  two  immiscible  fluids  as 
in  the  case  of  bubbles  or  drops  immersed  in  an  ambient  fluid. 

Whether  the  immersed  interface  represents  a  phase  interface  or  a  boundary  be- 
tween two  fluids  with  different  transport  properties,  the  numerical  simulation  of  these 
processes  has  to  deal  with  moving  interfaces.  The  numerical  investigation  of  moving 
interface  problems  has  posed  a  difficult  challenge  for  a  long  time,  especially  when  a 
moving  phase  change  interface  and  extremely  large  ratios  of  fluid  properties  across  the 
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Figure  1.5:  Illustration  of  multiple  immersed  interfaces  in  a  computational  domain 


interface  exist  simultaneously,  as  in  the  present  study  of  the  liquid- vapor  phase  change 
problem.  Despite  its  difficulties,  in  the  past  several  decades,  significant  advances  have 
been  made  in  numerical  techniques  for  obtaining  the  solutions  of  immersed,  moving 
interface  problems. 

In  this  section,  we  will  briefly  examine  the  available  numerical  approaches  for 
immersed,  moving  interface  problems  with  fluid  flow. 

Numerical  techniques  for  immersed,  moving  interface  problems  could  be  classi- 
fied into  three  broad  categories  (see  Shyy,  Udaykumar,  Rao  k  Smith,  1996):  body- 
fitted  grid  approach,  fixed  grid  approach  and  mixed  Eulerian-Lagrangian  approach. 

The  natural  approach  is  to  use  body-fitted  type  mesh  which  could  be  either  struc- 
tured or  unstructured  (e.g.,  finite  element  mesh).  Relatively  complex  two-dimensional 
free  boundary  problems  were  simulated  using  this  approach  (see  Ryskin  &  Leal,  1984; 
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Kang  &  Leal,  1987;  Dandy  &  Leal,  1989).  A  body-conforming  mesh  makes  the  bound- 
ary conditions  specification  on  the  interface  relatively  easy  and  accurate  because  the 
interfacial  conditions  are  applied  at  the  exact  location  of  the  interface  without  any 
smearing  operation.  The  sharp  discontinuity  of  material  properties  across  the  inter- 
face is  always  maintained.  However,  in  the  body-fitted  approach,  the  computational 
expense  is  an  issue  because  the  grid  has  to  be  modified  or  regenerated  continuously 
to  account  for  the  movement  of  the  interface.  The  body-fitted  approaches  encounter 
major  difficulty  in  dealing  with  the  topological  complexity  resulted  from  arbitrary 
interface  movement  with  desired  mesh  quality.  Besides,  the  governing  equations  are 
not  in  a  simple  form  for  body-fitted  mesh  which  causes  extra  difficulty  for  the  over- 
all solution  procedure.  Multiple  connectivity  makes  designing  an  efficient  algebraic 
solver  a  difficult  task. 

Another  completely  different  approach  is  to  use  the  fixed  Cartesian  grid.  In 
this  purely  Eulerian  approach,  the  exact  location  of  the  interface  is  not  explicitly 
known,  in  contrast,  the  interface  details  are  determined  from  solving  a  field  variable. 
Within  this  approach,  there  are  a  few  variants,  for  example,  the  Volume  of  Fluid 
(VOF)  method  (see  Hirt  &  Nichols,  1981;  Brackbill,  Kothe  &  Zemach,  1992;  Kothe 
&  Mjolsness,  1992;  Rider  &  Kothe,  1998;  Scardovelli  &  Zaleski,  1999)  ,  the  level-set 
method  (see  Osher  &  Sethian,  1988;  Sethian,  1990;  Sussman,  Smereka  &  Osher,  1994; 
Sussman,  Fatemi,  Smereka  &  Osher,  1998)  and  phase-field  method  (see  Caginalp, 
1984;  Kobayashi,  1993;  Antanovskii,  1995;  Anderson  &  McFadden,  1997;  Jacqmin, 
1999)  amongst  others. 

When  surface  tension  on  the  interface  can  not  be  ignored,  the  accurate  estimation 
of  the  interface  curvature  is  crucial  to  successfully  predicting  the  interface  evolution. 
Under  these  circumstances  where  explicit  interface  tracking  is  needed,  the  mixed 
Eulerian-Lagrangian  approach  seems  attractive  because  in  this  approach,  apart  from 
the  fixed  structured  grids  (Eulerian  part)  for  the  flow  field,  interfaces  are  explicitly 
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described  by  curves  or  surfaces  using  computational  geometry  tools  (Lagrangian  part) 
and  the  movement  of  those  curves  and  surfaces  are  determined  by  the  computational 
process. 

In  the  fixed  Cartesian  based  approach,  the  grid  generation  is  greatly  simplified 
and  the  governing  equations  retain  the  relatively  simple  form  in  Cartesian  coordi- 
nates. Because  the  grid  is  Cartesian,  many  efficient  solvers  are  easily  obtained  to 
reduce  the  computational  cost. 

In  the  mixed  approach,  a  major  issue  is  how  to  handle  the  interaction  between 
a  moving  interface  and  a  flow  field.  The  interface  would  cut  through  the  underlying 
Cartesian  grid  in  an  arbitrary  manner  which  resulted  in  fragmented  control  volumes 
(cut  cells)  around  the  interface.  Fragmented,  small  cut  cells  would  cause  numerical 
instability.  To  circumvent  this  difficulty,  the  so  called  immersed  boundary  technique 
has  been  widely  used  to  simulate  multi-fluid  flows  (see  Peskin,  1977;  Unverdi  k 
Tryggvason,  1992;  Goldstein,  Handler  k  Sirovich,  1995;  Juric  k  Tryggvason,  1996, 
1998;  Kan,  Udaykumar,  Shyy  k  Tran-Son-Tay,  1998).  In  this  technique,  the  interface 
is  not  kept  sharp  but  is  represented  on  the  fixed  underlying  mesh  by  a  thin  layer  of 
finite  thickness  around  the  interface.  A  scalar  indicator  function  is  then  constructed 
from  the  known  position  of  the  interface  to  smooth  out  the  jumps  in  the  material 
properties  across  this  thin  layer.  The  information  transfer  between  the  interface  and 
a  flow  field  is  done  by  introducing  additional  source  terms  in  the  governing  equations. 
This  technique  is  quite  robust,  however  its  fuzzy  interface  treatment  does  not  reflect 
closely  the  physical  reality  of  the  zero  thickness  interface,  and  thus  it  is  not  suitable 
for  accurately  resolving  any  boundary  layer  type  of  behavior  in  the  thermal  and  flow 
fields  when  convection  is  important. 

A  strictly  conservative  treatment  of  small  cut  cells  is  the  so-called  cut-cell  tech- 
nique. In  this  treatment,  typically,  small  cut  cells  are  absorbed  into  neighboring  cells 
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to  form  mostly  quadrilateral  cells  in  the  vicinity  of  the  interface.  Thereby,  the  inter- 
face is  treated  as  a  sharp  discontinuity  and  the  conservation  property  of  the  solution 
procedure  is  always  retained.  The  interface  effect  is  taken  into  account  by  the  flow 
solver  through  the  reshaped  control  volumes  near  the  interface. 

The  finite  volume  cut  cell  techniques  have  been  used  in  dealing  with  the  complex 
solid  boundaries  in  the  flow  field  in  conjunction  with  Cartesian  mesh  (see  Quirk, 
1994;  Udaykumar  &  Shyy,  1995;  Pember,  Bell,  Colella,  Crutchfield  &  Welcome,  1995; 
Udaykumar,  Kan,  Shyy  &  Tran-Son-Tay,  1997;  Ye,  Mittal,  Udaykumar  &  Shyy, 
1999).  In  those  applications,  the  relative  locations  between  solid  boundaries  and 
grid  cells  are  fixed  throughout  the  solution  process.  There  is  no  need  to  determine 
the  boundary  shape  in  the  solution  process.  There  is  no  interaction  between  rigid 
boundaries  and  flow  field  which  makes  it  much  easier  to  implement  the  cut  cell 
technique. 

In  the  applications  where  free,  moving  interfaces  with  phase  change,  e.g.,  in- 
terfaces separating  bubbles  or  drops  from  the  surrounding  fluids,  are  involved.  It 
is  much  more  challenging  to  implement  the  cut  cell  technique  because  the  interface 
shape  is  not  fixed  any  more,  rather  depending  on  the  dynamic  conditions  on  the 
interface.  In  other  words,  the  instantaneous  interface  shape  must  be  determined  as 
part  of  overall  solution  process  based  on  the  interaction  of  the  interface  with  the  flow 
fields  in  both  sides  of  the  interface.  Owing  to  this  requirement,  no  reports  has  been 
found  in  the  literature  on  the  application  of  finite  volume  cut  cell  techniques  to  free 
interfaces  like  those  of  bubbles  or  drops.  In  the  present  research,  a  fixed  grid,  sharp 
interface  moving  boundary  method  is  developed  and  applied  to  simulate  the  moving 
bubble  interfaces. 

To  summarize,  pure  body-fitted  grid  approaches  are  seldom  used  nowadays,  they 
are  more  often  combined  with  fixed  grid  methods  to  improve  the  efficiency.  In  such 
configurations,  most  of  the  domain  is  covered  by  the  fixed  grid,  body-fitted  grid  is 
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only  used  around  the  curved  interface  to  capture  the  movement  of  the  interface.  A 
communication  mechanism  has  to  be  designed  between  global  fixed  grid  and  local 
body-fitted  grid. 

Fixed  grid  methods  are  mainly  aimed  at  industrial  simulations.  In  these  meth- 
ods, normally  a  homogeneous  system  is  applied  to  the  entire  domain.  As  a  result,  a 
single  set  of  governing  equations  are  solved  in  a  domain  including  multiphase  with 
different  properties.  The  implementation  of  these  methods  is  relatively  easy.  Despite 
its  simplicity,  a  sharp  discontinuity  for  the  transport  properties  at  the  phase  front 
should  be  preserved  in  conformity  with  the  physical  reality.  In  fixed  grid  methods, 
the  interface  is  not  explicitly  known  and  has  a  finite  thickness,  or  sometimes,  al- 
though the  interface  itself  can  be  implicitly  captured  by  some  field  function  to  be 
sharp,  large  property  jumps  across  the  two  phases  still  have  to  be  smoothed  over  a 
few  grid  points  around  the  interface  to  overcome  the  numerical  stiffness  caused  by 
the  property  jumps.  In  order  to  maintain  a  sharp  front,  the  fuzzy  interface  should  be 
restricted  to  within  a  few  grid  cells:  the  less  the  grid  cell,  the  better  the  resolution 
of  the  interface. 

Our  research  is  intended  to  capture  the  local  physics  near  the  phase  interface  for 
liquid-vapor  phase  change  process.  We  developed  a  numerical  method  that  is  based 
on  fixed  grid  with  explicitly  tracked  interfaces.  At  the  same  time,  the  explicitly 
tracked  interface  is  a  sharp  discontinuity  of  fluid  properties.  The  current  method  is 
a  mixed  Eulerian-Lagrangian  approach.  However,  the  major  difference  between  this 
method  and  other  mixed  Eulerian-Lagrangian  methods  is  the  interface  treatment. 
The  sharp  interface  is  maintained  in  the  present  method. 

In  the  present  method,  first,  an  explicit  interface  representing  and  tracking  pro- 
cedure is  developed  to  locate  the  interface  exactly,  and  a  underlying  fixed  grid  is 
employed  to  cover  the  entire  domain  for  field  solution.  Secondly,  the  interface  is  kept 
sharp  in  terms  of  property  jumps.  Thirdly,  two  set  of  governing  equations  are  solved 
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in  two  phases,  respectively  in  conjunction  with  corresponding  transport  properties. 
As  a  result,  the  fixed  grid  over  the  entire  domain  is  virtually  divided  into  two  sub- 
domains  using  cut  cell  techniques.  The  coupling  of  the  two  phases  is  treated  with 
the  implicit  updating  of  the  interface  location  under  the  constraints  of  interfacial 
conditions. 

1.3     Research  Objectives 
1.3.1     Research  Objectives 

The  interests  in  the  vapor  bubble  phenomena  has  stimulated  continuous,  inten- 
sive efforts  from  numerous  investigators  on  this  subject  for  half  a  century,  as  described 
in  the  previous  sections.  And  liquid-vapor  phase  change  (boiling  or  condensation) 
process  is  one  of  the  most  difficult  heat  transfer  modes  to  investigate.  Because  the 
vapor  bubble  behavior  is  described  by  a  set  of  coupled,  nonlinear  partial  differential 
equations  and  because  there  are  not  mathematical  tools  to  solve  those  equations  an- 
alytically, the  majority  of  the  theoretical  works  on  the  vapor  bubble  dynamics  in  the 
literature  have  to  resort  to  more  or  less  simplified  theoretical  models  so  as  to  make 
the  problem  tractable.  Although  those  simplified  theoretical  models  sometimes  do 
yield  quite  good  results  in  comparison  with  experimental  data,  it  is  certainly  desired 
and  has  been  ceaseless  attempted  to  have  the  theoretical  model  as  close  as  possible 
to  that  set  of  coupled,  nonlinear  partial  differential  equations.  The  present  research 
thus  serves  to  bring  that  closeness  to  a  level  which  has  not  been  achieved  before, 
however,  the  gap  between  the  theoretical  model  and  the  ultimate  description  is  not 
completed  removed. 

On  the  other  hand,  owing  to  experimental  difficulties,  only  a  few  investigations 
have  focused  on  the  micro-phenomena  at  the  phase  change  interface.  In  order  to 
achieve  a  comprehensive  understanding  of  the  physics  of  boiling  and  condensation, 
an  effort  is  needed  in  the  development  of  a  direct  numerical  simulation  technique  that 
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can  be  used  to  study  the  microscopic  transport  mechanisms  at  the  moving  liquid- 
vapor  phase  change  interface. 

The  objective  of  the  present  research  is  thereby  to  develop  a  numerical  method  to 
simulate  the  liquid-vapor  phase  change  transport  mechanisms  of  a  translating  vapor 
bubble  in  a  superheated  or  subcooled  liquid.  Based  on  the  numerical  simulation 
results,  insights  are  offered  into  the  physical  and  numerical  characteristics  of  the 
bubble  dynamics.  An  assessment  of  some  conclusions  or  correlations  proposed  in  the 
literature  regarding  the  bubble  growth  or  collapse  dynamics  is  also  presented. 

The  development  of  the  present  fixed  grid,  sharp  interface  moving  boundary 
method  consists  of  the  implementations  of  the  following  major  components: 

-  Interface  representation  algorithms  using  B-spline  curve  fitting, 

-  Fairing  algorithm  for  continuous  curvature  calculations, 

-  Intersection  of  the  curved  interface  with  grid  lines, 

-  Assemble  and  reshaping  of  the  cut  cells, 

-  Discretization  of  governing  equations  on  regular  and  cut  cells, 

-  Fractional  step  procedure  for  solving  the  coupled  governing  equations, 

-  Imposition  of  the  boundary  conditions  for  velocity,  pressure  and  temperature 
at  the  interface, 

-  Iterative  algorithm  for  updating  the  interface  location  based  on  normal  stress 
balance  conditions  at  the  interface, 

-  Algorithm  for  the  interface  evolution  owing  to  phase  change  effect. 

This  challenging  task  is  accomplished  owing  to  the  latest  development  of  numeri- 
cal techniques  for  solving  time-dependent  free  boundary  problems  in  fluid  mechanics. 
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In  recent  years,  it  has  become  possible  to  solve  coupled,  nonlinear  mass,  momen- 
tum and  energy  conservation  equations  (partial  differential  equations)  numerically 
because  advanced  numerical  techniques  and  fast  computer  hardware  were  available. 
For  problems  with  even  more  physical  complexity  such  as  those  involving  multi-phase 
fluid  flow,  heat  and  mass  transfer  and/or  moving  interfaces  etc.,  the  reported  simu- 
lation works  on  Navier-Stokes  equations  start  to  appear  in  the  literature  only  in  the 
last  ten  years  yet  there  are  many  numerically  relevant  issues  which  are  not  resolved. 

In  our  case,  the  method  for  solving  the  free  boundary  problems  is  developed  and 
large  scale  computations  are  carried  out  to  solve  coupled  Navier-Stokes  and  energy 
equations  in  liquid  and  vapor  phase  respectively  at  the  same  time.  On  the  phase 
transition  interface,  the  mass,  momentum  and  energy  conservation  principles  are 
matched  from  the  liquid  and  vapor  sides  of  the  interface. 

The  method  developed  in  the  current  research  is  not  restricted  to  only  liquid- 
vapor  phase  change  problems,  it  is  generally  applicable  to  other  phase  change  pro- 
cesses such  as  solid-liquid  phase  change  of  melting  or  solidification  because  the 
method  is  developed  to  handle  property  discontinuity  and  moving  interfaces  which 
are  common  to  the  liquid-vapor  and  solid-liquid  phase  change. 

1.3.2     Scales  and  Dimensionless  Parameters 

Several  major  non-dimensional  parameters  including  Reynolds,  Weber,  Jakob 
and  Fourier  numbers  are  varied  to  gain  a  better  understanding  of  physical  and  nu- 
merical aspects  of  the  solution  procedure.  To  determine  the  above  dimensionless  pa- 
rameters, characteristic  scales  are  chosen  based  on  the  problem  features.  The  length 
and  velocity  scales  are  determined  as  well  as  the  maximum  temperature  difference, 
AT,  in  the  system.  Naturally,  we  choose  the  initial  bubble  diameter,  D0,  as  the 
length  scale  L.  The  velocity  scale  can  be  defined  based  on  the  diffusion  mechanism: 
OLijL  for  stationary  bubble  growth  cases,  buoyancy  effect:   (gL)ll2  for  rising  bubble 
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Table  1.3:  Saturation  properties  of  water  at  Psat  =  101.3  kPa,  Tsat  —  373.15 /f 


Saturation  Properties 

Liquid 

Vapor 

P  (kg/m3) 

958.3 

0.597 

fi  (N  s/m2) 

277.53  x  10"6 

12.55  x  10"6 

cp{J/kgK) 

4,220 

2,030 

k{W/mK) 

0.679 

0.025 

Pr 

1.72 

1.02 

A  (J/ kg) 

2,256,700 

a  (N/m) 

0.0589 

Table  1.4:  Saturation  properties  of  Ethanol  at  Psat  =  101.3  kPa,  Tsat  —  351.45  K 


Saturation  Properties 

Liquid 

Vapor 

P  (kg/m3) 

757.0 

1.435 

fx  (N  s/m2) 

428.7  x  10"6 

10.4  x  10~6 

cp(J/kgK) 

3000 

1830 

k{W/mK) 

0.1536 

0.0199 

Pr 

8.37 

0.96 

A (J /kg) 

963, 

000 

a  (N/m) 

0.0177 

with  phase  change  cases,  or  the  bubble  terminal  velocity  U  for  bubble  rising  with  no 
phase  change.  The  choice  in  each  case  will  be  specified  individually. 

The  rest  of  the  information  required  to  specify  the  dimensionless  parameters  are 
the  thermo-physical  properties  of  the  liquid.  An  example  of  working  fluids  is  water 
whose  saturation  properties  at  the  atmospheric  pressure  are  listed  in  Table  1.3.  Then 
the  typical  magnitude  of  those  dimensionless  property  ratios  of  the  liquid  and  vapor 
phases  of  water  are  pi/pv  =  1605.2,  m/fa  =  22.114,  kt/kv  =  27.16,  cPl/cPv  =  2.0788. 

The  range  of  Reynolds  number  we  consider  in  this  research  is  Re  =  1  ^  1000. 
The  bubble  diameter  is  0.5  mm  >-.  1.5  mm,  and  the  maximum  temperature  difference 
AT  =  5°C. 

Besides  water,  we  also  consider  another  working  fluid,  namely,  Ethanol  (CH3CH2OH) 
whose  thermo-physical  properties  are  listed  in  Table  1.4. 


CHAPTER  2 
NUMERICAL  METHOD 

In  this  chapter,  the  fixed-grid,  sharp  interface  moving  boundary  method  is  de- 
scribed. In  this  method,  a  underlying  fixed  Cartesian  grid  serves  as  the  Eulerian 
framework  of  the  algorithm.  Within  this  Eulerian  framework,  separate  marker  points 
which  are  connected  by  piece-wise  polynomials  are  adopted  to  represent  the  phase 
interface.  Those  regions  separated  by  the  phase  interface  represent  different  phases. 
The  deformation  and  movement  of  the  interface  is  realized  through  the  translation  of 
these  markers  over  the  underlying  stationary  grid,  which  constitutes  the  Lagrangian 
part  of  the  algorithm.  In  each  phase  region,  the  fractional  step  method  is  employed 
to  solve  the  coupled  governing  equations  of  either  liquid  or  vapor  phase. 

2.1     Fixed-Grid,  Sharp  Interface  Method 

The  development  of  this  numerical  method  consists  of  three  major  tasks: 

•  Implementation  of  curve  fitting  and  curvature  calculation  algorithms. 

•  Development  of  Cartesian  grid  techniques  for  solving  fluid  flows  over  fixed, 
immersed  interfaces  with  complex  geometries. 

•  Development  of  moving  interface  algorithms,  such  as  the  interface  updating 
algorithm. 

The  content  of  this  chapter  is  organized,  as  shown  in  Table  2.1,  so  as  to  describe 
these  three  major  components  of  the  numerical  method  in  a  logic  manner.  In  the 
section  2.2  of  this  chapter,  the  implementation  of  interface  representing,  tracking  and 
curvature  calculation  algorithms  are  detailed.    The  C2  cubic  B-spline  curve  fitting 
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Table  2.1:  Organization  of  Chapter  2 


Interface  representing  and  tracking  algorithms 

Section  2.2 

C2  cubic  B-spline  algorithm 

Fairing  algorithm 

Validation  of  geometric  algorithms 

Section  2.2.1 
Section  2.2.2 
Section  2.2.3 

Cartesian  grid  method  for  fixed,  immersed  interfaces 

Section  2.3 

Fractional  step  procedure  for  solving  coupled  governing  equations 
Discretization  of  governing  equations  on  regular  and  cut  cells 
Imposition  of  boundary  conditions  at  the  interface 
Algebraic  equation  solver 

Section  2.3.1 
Section  2.3.2 
Section  2.3.3 
Section  2.3.4 

Moving  interface  algorithms 

Section  2.4 

Reinitialize  the  flow  field  due  to  phase  change 
Inclusion  of  interfacial  conditions 
Iterative  Process  for  Determining  the  Interface  Shape 
Determining  Interface  Shape  with  Phase  Change 
Discontinuity  of  Material  Properties 

Section  2.4.1 
Section  2.4.2 
Section  2.4.3 
Section  2.4.4 
Section  2.4.5 

is  employed  in  conjunction  with  fairing  algorithm.  In  the  section  2.3,  numerical 
procedures  related  to  solving  fluid  flows  over  fixed,  immersed  interfaces  on  Cartesian 
grid  are  described.  In  the  section  2.4,  the  numerical  techniques  for  moving  interface 
computations  are  presented. 

The  overall  solution  procedure  for  numerical  simulation  of  moving  interface  prob- 
lems is  illustrated  in  the  flow  chart  of  the  code  in  Figure  2.1  in  which  the  major 
function  modules  interface  tracking  procedure  is  detailed  in  Figure  2.2,  fractional  step 
process  is  detailed  in  Figure  2.3,  moving  and  updating  interface  is  detailed  in  Figure 

2.4.  In  the  rest  of  the  chapter,  each  section  is  dedicated  to  certain  numerical  issues 
appeared  in  these  flow  charts  as  also  shown  in  Table  2.1. 

Currently,  the  source  code  of  this  solver  has  nearly  40,000  lines  of  Fortran90 
statement  implemented  on  Compaq  Tru64  Unix  work  stations. 

2.2     Interface  Representation  Algorithms 

A  representative  schematic  of  the  present  fixed  grid  method  is  depicted  in  Figure 

2.5.  The  grid  is  Cartesian  and  does  not  conform  to  the  body  surface,  and  the  interface 
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(   Begin  J 
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Input  computational  parameters 
Input  initial  interface  shape 


I 


Interface  tracking  procedure  (Section  2.2) 


Time  step  n=0 


\ 


Time  step  loop 
Time  step  n+1 


•  Moving  and  updating  interface  (Section  2.4) 

•  Interface  tracking  procedure  (Section  2.2) 

•  Cartesian  Grid  solver  (Section  2.3) 


Until  time  step  >  specified  number 


X 

L_EndJ 


Figure  2.1:  The  flow  chart  showing  the  overall  solution  procedure  of  moving  interface 
algorithm. 
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Redistribute  the  marker  points  to  have  equal  spacing 


Parameterization  of  the  interface  through  marker  points 


Find  the  intersecitons  of  the  interface  with  grid  lines 


Identify  the  phase  which  each  cell  lies  in  and  flag  cells 


Creat  and  assemble  the  cut  cells 


Assign  the  material  properties 


Figure  2.2:  The  flow  chart  showing  the  interface  tracking  procedure  as  detailed  in 
Section  2.2. 
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Initialize  the  flow  field 


Solve  for  intermediate  flow  field  u*  (Eqn.  2.12) 

•  Imposition  of  boundary  conditions  (Section  2.3.3) 

•  Calculate  face  fluxes  using  interpolation  procedure  (Section  2.3.2) 


I 


Obtain  face  center  velocity  U* 


I 


Solve  for  pressure  Poisson  equation  (Eqn.  2.15) 

•  Imposition  of  boundary  conditions  (Section  2.3.3) 

•  Calculate  face  fluxes  using  interpolation  procedure  (Section  2.3.2) 


n+l 


Correct  cell  center  velocity  field  u*—>  u 
Correct  cell  center  velocity  field  U*—>  U"*1 
(Eqn.  2.16) 


I 


Solve  for  energy  equation  T" 


+i 


Figure  2.3:  The  flow  chart  showing  the  fractional  step  procedure  as  detailed  in  Section 
2.3. 
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Iteration  k=0 

1 

Interface  update  iteration 
Iteration  k=k+l 

1 

Move  interface  to  new  location  (Section  2.4.3) 

1 

Modify  the  interface  shape  to  satisfy  mass 
conservation  in  each  phase  (Section  2.4.3) 

t 

Interface  representation  procedure  (Section  2.2) 

t 

Reinitialize  the  flow  field  due  to  phase  change  (Section  2.4.1} 
Fractional  step  procedure  (Section  2.3.1) 

1 

Terminated  if  all  conditions  are  satisfied 

Fig 
det 

lire  2.4:  The  flow  chart  showing  the  interface  moving  and  updating  procedure 
ailed  in  Section  2.4. 
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Figure  2.5:  Schematic  of  Cartesian  grid  and  interfacial  marker  points. 


is  defined  explicitly  by  geometric  curves  in  the  computational  domain.  Basic  elements 
involved  in  defining  the  interface  are  interface  marker  points  and  interfacial  curves 
connecting  the  markers.  The  markers  define  the  terminal  points  of  the  interfacial 
curves.  Given  a  set  of  markers,  finding  a  curve  which  fits  all  markers  is  a  geometric 
interpolation  process.  Depending  on  geometric  conditions  imposed  at  the  marker 
points,  there  are  various  ways  to  define  numerically  the  interface  characteristics. 

In  the  present  work,  the  C2  cubic  interpolatory  B-spline  curve  is  employed  to 
fit  a  set  of  marker  points.  The  B-spline  curve  refers  to  a  set  of  Bezier  curves  glued 
together  at  the  marker  points  to  form  a  composite  curve  representing  the  interface. 
Each  piecewise  Bezier  curve  is  actually  a  polynomial  curve  expressed  mathematically 
for  convenience  in  special  form:  Bernstein  polynomials.  There  are  efficient  algorithms 
to  do  curve  manipulation  and  geometric  calculations.  CT  denotes  the  smoothness  con- 
ditions at  the  junction  points  of  piecewise  Bezier  curves,  requiring  that  a  composite 
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curve  is  r  times  continuously  differentiable  at  the  junction  points.  C2  thereby  means 
the  composite  curve  has  continuous  second  derivatives  at  any  marker  point.  In  order 
to  achieve  this  level  of  smoothness,  the  individual  Bezier  curve  must  be  at  least  degree 
3  (cubic)  polynomial  to  have  continuous  second  derivatives,  and  second  derivatives 
computed  from  both  sides  of  the  junction  point  have  to  be  equal.  This  mathematical 
requirement  of  the  two  polynomials  yields  a  smooth  global  curve  to  represent  the 
entire  interface. 

Once  the  B-spline  fitting  curve  of  the  interface  is  constructed,  the  geometric 
informations  such  as  location  and  curvature  of  any  point  along  the  entire  interface 
can  be  easily  obtained. 

The  translating,  deforming,  expanding  or  shrinking  of  the  interfaces  are  realized 
through  the  motion  of  each  individual  marker  point  which  in  turn  is  determined  from 
the  flow  quantities  on  underlying  fixed  grid  using,  e.g.,  normal  stress  balance  condi- 
tions. With  the  movement  of  the  markers,  the  instantaneous  B-spline  representation 
of  the  interface  is  reconstructed  accordingly  to  keep  track  of  the  interface.  For  a 
interface  in  motion,  the  curvature  calculated  based  on  the  continuously  constructed 
B-spline  may  exhibit  numerical  oscillations  (see  Farin,  1997).  To  extract  the  correct 
curvature  values  along  the  interface,  a  so-called  fairing  algorithm  is  adopted.  As  will 
be  presented,  the  combination  of  the  cubic  B-spline  and  the  fairing  algorithm  results 
in  a  robust  and  accurate  method  for  tracking  highly  distorted  interfaces  in  terms  of 
location  as  well  as  curvature. 

2.2.1     C2  Cubic  B-spline  Algorithm 

The  interface  representation  using  B-spline  is  based  on  the  marker  locations  com- 
puted at  every  time  instant.  The  interfacial  marker  points  are  indexed  sequentially 
and  can  represent  any  number  of  open  or  closed  interfaces  based  on  the  assigned 
connectivity  between  them.  There  is  no  restriction  on  the  number  of  interfaces  or 
on  the  end  conditions  that  can  be  imposed  on  these  interface  curves.   As  shown  in 
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X,<s) 


Phase  0 


Phase  1 


Figure  2.6:  Illustration  of  immersed  interfaces,  marker  points  and  normal  convention. 


Figure  2.6,  the  interfaces  can  be  open  or  closed,  stationary  or  moving  and  can  enclose 
different  phases  or  materials  with  different  properties.  The  locations  of  the  marker 
points  are  defined  by  coordinates  X(s)  which  is  parameterized  in  terms  of  arc  length 
s  and  a  distance  ratio  based  on  the  grid  spacing.  In  our  investigation,  the  marker 
spacing  is  initially  assigned  to  be  the  same  as  the  grid  spacing,  h.  In  the  course  of 
computation,  markers  are  redistributed  after  each  time  step  according  to  the  initial 
criterion. 

The  convention  adopted  by  the  present  algorithm  in  indexing  the  marker  points 
is  that,  by  choosing  the  normal  of  the  interface  to  point  from  phase  1  to  phase  0,  the 
phase  1  always  lies  to  the  right  as  one  traverses  the  interface  along  the  sequence  of 
the  marker  points. 

In  order  to  obtain  geometric  information,  we  need  to  construct  a  smooth  curve 
fitting  this  set  of  marker  points.  In  the  computer-aided  geometric  design  applications, 
the  B-spline  offers  the  greatest  flexibility  and  effectiveness  to  serve  this  purpose.  The 
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B-spline  curves  are  composed  of  piece- wise  Bezier  curves  between  any  two  consecutive 
markers.  When  we  connect  each  Bezier  curve  piece  by  piece  to  form  the  global  B- 
spline  curve,  the  original  data  points  (here  marker  points)  are  junction  points  of  two 
consecutive  Bezier  curves.  Now  the  question  arises:  how  is  the  smoothness  of  this 
overall  B-spline  curve?  This  question  is  answered  by  the  definition  of  Cn  continuity 
condition  of  the  B-spline  curves.  The  Cn  continuity  means  that  at  any  junction  point, 
the  nth  derivative  computed  from  the  left  Bezier  curve  equals  to  that  computed  from 
right  Bezier  curve.  The  most  popular  B-spline  curve  is  the  C2  spline.  To  have  C2 
continuity,  each  Bezier  curve  must  be  at  least  degree  three;  therefore  the  cubic  Bezier 
curves  can  be  used  as  piece-wise  polynomial  in  this  context.  The  curve  fitting  scheme 
we  adopted  in  this  research  is  C2  cubic  B-spline  curve. 

In  the  following,  we  describe  the  construction  of  C2  cubic  B-spline  curves  and 
so-called  'fairing'  algorithm  which  is  used  to  removed  numerical  noise  to  recover  the 
accurate  information  of  the  curvature  of  the  interface. 

Based  on  the  spline  theory  (see  Ahlberg  et  ai,  1967;  Farin,  1997),  to  construct 
a  C2  cubic  B-spline  curve  interpolating  a  set  of  data  points  Xo,-..,Xl  with  corre- 
sponding parameter  values  (or  knots)  u0,...uL,  the  vertices  d{  of  the  B-spline  control 
polygon  have  to  be  determined  first,  as  shown  in  Figure  2.7.  In  this  example,  five 
points  x0,  ■  ■  ■  ,x5  are  assigned  initially.  The  corresponding  knot  sequence  u0,  ■  ■  ■  ,u5 
is  chosen  as  the  chord  length  at  each  point.  The  example  here  illustrates  a  closed 
(periodic)  curve,  so  x5  =  x0.  The  relationship  between  the  data  points  Xi  and  the 
control  vertices  di  is  (see  Farin,  1997) 

(Z\i_i  +  Ai)  Xi  =  Oidi-i  +  Pidi  +  7idi+1  (2.1) 
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Figure  2.7:   Cubic  interpolatory  spline  curve  of  five  points  along  with  its  B-spline 
control  polygon,  Bezier  control  polygon. 
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where  we  have  Ai  =  Aui  =  itj+i  —  Ui  and 

(A)2 


0L{ 


Ai-2  +  A-i  +  A 


AjjAj-2  +  A-i)    {   Aj^jAj  +  A+i) 
Pi        Ai_2  +  A-i  +  A     A-i  +  A  +  A+i 
(A-i)2 


li   = 


At-i  +  A  +  Ai+l 
Using  the  periodic  condition,  we  obtain  a  linear  system  of  the  form 

a0 


(2.2) 


A)      7o 
ai      ft    7i 
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OtL-l      $L-1 
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ro 

di 

ri 

dL-2 

TL-2 

dL-i 

.  rL~l   . 

(2.3) 


where  the  right-hand  sides  are  of  the  form 

n  =  {A-i  +  a^  Xi 

The  equation  Eq.  (2.3)  can  be  solved  using  the  procedure  described  in,  e.g., 
Ahlberg  et  al.  (1967),  P.  15,  yielding  the  vertices  d0,  •  •  •  ,  d5  as  shown  in  Figure  2.7. 

Then  the  vertices  of  control  polygon  for  piece-wise  cubic  Bezier  curves  can  be 
obtained  based  on  this  B-spline  control  polygon,  namely,  using  the  parameter  values 
to  determine  two  vertices  on  each  leg  of  the  B-spline  control  polygon  as  shown  in 
Figure  2.7.  Each  of  the  piece-wise  Bezier  curve  is  defined  by  four  control  vertices, 
including  the  two  marker  points,  denoted  by  solid  circles,  and,  between  them,  two 
vertices  to  be  decided  from  the  B-spline  control  polygon.  Therefore  all  four  Bezier 
curve  control  vertices  are 


&3i 

hi-2 

&3J-1 


A-\  +  Ai 


di_i  H -r—  di 


id-  + 


A-2  +  A-i 


dt 


(2.4) 
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where  A  —  A^2  +  ^i-i  +  A{.  The  index  i  denotes  the  ith  vertice  of  B-spline  control 
polygon,  and  the  total  number  of  Bezier  control  vertices  is  3z  corresponding  to  vertice 
i  of  B-spline  control  polygon  for  a  closed  curve. 

Once  we  have  the  control  vertices  for  each  local  Bezier  curve,  by  defining  a  local 
parameter  0  <  t  <  1  for  the  interval  [M{,Ui+i]  as  t  =  J*"H!Li  we  can  express  the 
piecewise  cubic  Bezier  curve  in  the  form  of 

3 

6(t)  =  5>fl?(t)  (2-5) 

i=0 

where  Bernstein  polynomial  Bf{t)  is  defined  by 

Bf (*)  -  (.J  )  **(1  -  *)*-*  (2-6) 

with  binomial  coefficient 

3  N  m  {  a^feji  »/  °  < » < 3 

i  )      \  0  e/se 

Using  Eq.  (2.5),  the  geometric  information  such  as  normal  and  curvature  etc. 
can  be  evaluated  using  analytical  formulas  as 

Vt 


nx    = 


ny 


K  "    W  +  yl?'2  {    ] 

For  axisymmetric  geometries,  the  total  curvature  is  the  sum  of  «  in  Eq.  (2.7) 
and  the  other  principle  curvature,  i.e.,  xt/y{x}  +  yt2)3^2. 

2.2.2     Fairing  Algorithm 

From  the  numerical  simulation  point  of  view,  the  smoothness  of  the  interfacial 
curve  and  the  smoothness  of  curvatures  are  two  most  important  characteristics  of 
the  interface.  Using  the  B-spline  fitting,  while  a  curve  looks  apparently  smooth,  the 
curvatures  obtained  by  Eq.    (2.7)  can  be  contaminated  by  numerical  noises.    The 
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Figure  2.8:  A  sample  interfacial  curve  from  one  simulation  case 


curvature  formula  involves  the  second  derivatives  as  well  as  nonlinear  products  of 
the  first  derivatives,  and  is  prone  to  suffer  from  numerical  noises  even  if  the  interface 
shape  and  slope  are  found  to  be  smooth.  This  usually  occurs  when  the  interface 
is  evolving  in  the  course  of  computation.  As  an  example  to  illustrate  the  problem, 
Figure  2.8  shows  an  interface  curve  taken  from  a  simulation  conducted  in  this  study. 
The  corresponding  curvature  plot  based  on  B-spline  fitting  process  is  shown  in  Figure 
2.9.  There  are  substantial  oscillations  in  the  curvature  profile,  indicating  that  noises 
associated  with  the  numerical  procedures  are  substantial.  Such  phenomena  are  well 
known  in  computer-aided  geometric  design  (see  Farin,  1997). 

To  treat  this  difficulty,  the  so-called  curve  fairing  (smoothing)  algorithms  have 
been  developed  in  the  literature  (see  Sapidis  &  Farin,  1990;  Farin,  1997).  One  such 
algorithm,  in  the  C2  cubic  spline  case,  is  to  make  local  spline  curve  segment  three 
times  differentiable  which  is  one  order  higher  than  that  required  by  the  original 
cubic  spline.    The  basic  idea  is  to  adjust  the  vertex  locations  to  make  the  local 
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Figure  2.9:  Curvature  plot  obtained  from  B-spline  fitting  for  the  interfacial  curve  in 
Figure  2.8. 


(2.8) 


curve  segments  around  them  become  three  times  differentiable.  Since  the  curvature 
involves  only  first  and  second  derivatives,  this  extra  differentiable  requirement  will 
ensure  that  the  curvature  be  differentiable,  and  prevents  discontinuity  in  the  slope 
of  the  curvature,  as  shown  in  Figure  2.10.  The  formula  for  obtaining  new  B-spline 
control  vertex  dj  are  briefly  listed  as 

Uj   — 

uj+2  ~  Uj-2 

where  the  auxiliary  points  lj,  r±  are  given  by 

__    (uj+i  -  «,-_3)dj-i  -  (uj+i  -  uAdj-a 

ij    — 

Uj  —  Uj-3 
r        =       (Uj+3  ~  uj-l)dj  +  l  -  {Uj  -  Uj-i)dj+2 

uj+3  -  Uj 

The  geometry  interpretation  is  illustrated  in  Figure  2.10  which  is  from  page  365  of 
Farin  (1997).  The  detailed  discussions  can  be  found  in  Sapidis  &  Farin  (1990);  Farin 
(1997). 


(2.9) 
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Figure  2.10:  Knot  removal:  if  dj  is  moved  to  dj,  the  new  curve,  denoted  by  dotted 
line,  is  three  times  differentiable  at  Uj  [page  365,  Farin  (1997)]. 


In  practice,  typically  one  fairing  operation  is  not  sufficient,  and  the  fairing  pro- 
cedure will  be  repeated  multiple  times.  The  criterion  adopted  in  the  present  study  for 
determining  the  convergence  of  the  curvature  profile  is  that  the  maximum  difference 
of  the  curvature  values  among  all  marker  points  between  two  consecutive  fairing  op- 
erations must  be  less  than  10~3.  For  the  curve  in  Figure  2.8,  the  resulting  curvature 
plot  after  100  iterations  of  fairing  operation  is  shown  in  Figure  2.11.  The  maximum 
correction  in  the  curvature  is  within  10~3. 

It  is  noted  that  the  fairing  algorithm  is  a  geometric  operation.  It  obtains  correct 
curvatures  by  modifying  the  geometry  of  the  interface  itself,  not  by  manipulating 
the  formulas  for  computing  the  curvature  of  a  given  geometry.  On  the  other  hand, 
attention  needs  to  be  paid  to  ensure  that  the  interfacial  marker  locations,  as  well  as 
the  volume/surface  enclosed  by  the  interface,  be  satisfactorily  preserved.  A  critical 
criterion  of  developing  a  satisfactory  fairing  algorithm  is  that  with  arbitrary  number 
of  fairing  iterations,  the  geometric  information  can  be  maintained  at  an  asymptot- 
ically constant  state  without  being  continuously  smeared.  The  corrections  in  the 
interfacial  locations  are  usually  very  small,  so  the  corrected  marker  locations  and  the 
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Figure  2.11:  Curvature  plot  obtained  before  and  after  fairing  operation  for  the  curve 
in  Figure  2.8. 


resulted  interface  shape  will  not  have  perceptible  changes  after  fairing  operation,  as 
will  be  demonstrated  in  the  testing  of  geometric  algorithms. 

2.2.3     Validation  of  Geometric  Algorithms 

The  geometric  algorithms  for  interface  representing  and  tracking  described  in 
the  section  2.2.1  and  2.2.2  are  capable  of  handling  closed  and  open  interfaces.  They 
are  validated  here  for  accuracy  on  some  cases  used  in  Shyy  (1994). 

First  we  examine  the  expansion  of  a  circle  on  a  162  x  162  grid  as  in  Figure  2.12. 
Initially,  the  circle  has  a  radius  of  1,  and  expands  at  a  constant  speed.  The  curvature 
of  the  circle  at  any  given  time  is  a  constant.  Figure  2.12(a)  shows  the  interfacial 
shapes  at  equal  intervals  of  time.  Figure  2.12(b)  shows  the  curvature  computed 
by  the  aforementioned  formulas,  at  the  corresponding  time  instants  as  in  Figure 
2.12(a).  As  shown  there,  the  interfacial  curvature  is  a  constant  while,  by  regulating 
the  spacing  between  two  neighboring  markers,  the  number  of  markers  increases  as  the 
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circle  expands.  The  results  shown  in  Figure  2.12  are  obtained  from  B-spline  fitting 
algorithm  only.  No  fairing  operation  is  needed  for  this  case;  by  applying  fairing,  no 
impact  is  observed,  either. 

As  another  example  for  variable  curvatures,  Figure  2.13  shows  an  initial  interfa- 
cial  curve  shape  expressed  in  the  form  of  F,  =  2.0+ (1  —  cos(2irXi/10)).  The  curvature 
plot  is  found  to  be  smooth  using  B-spline  fitting  for  the  initial  shape  as  shown  in 
Figure  2.13(b).  Given  the  normal  velocity  of  each  marker  point  as  Vni  =  Yi  —  Yi,  the 
interface  moves  as  depicted  in  Figure  2.16(a).  Figure  2.14  shows  the  interface  shape 
and  corresponding  curvatures  at  the  early  stage  of  the  interface  evolution.  If  the 
fairing  is  not  applied,  the  curvature  obtained  from  B-spin  representation,  in  Figure 
2.14(b),  exhibits  oscillations  near  both  ends  of  the  curve  after  the  interface  evolves 
to  a  new  shape.  After  10  fairing  operations,  the  shape  and  curvatures  at  the  same 
time  instant  are  shown  in  Figure  2.15.  The  curvature  profile  in  Figure  2.15(b)  is  free 
from  the  artificial  spikes  at  the  ends. 

The  curvature  plot  in  Figure  2.15(b)  becomes  smooth  and  the  shape  in  Figure 
2.15(a)  is  well  preserved  and  can't  be  distinguished  from  that  in  Figure  2.14(a)  as  if 
the  shape  is  identical  to  previous  one.  This  is  further  confirmed  by  plots  in  Figure 
2.18.  Figure  2.18  shows  the  corrections  on  coordinates  of  markers  by  the  fairing 
operation  at  various  time  instants  from  the  beginning  to  the  end  of  the  development 
of  the  interface.  All  corrections  of  coordinates  are  small,  which  is  the  reason  why 
the  two  shapes  before  and  after  the  fairing  operations  are  virtually  identical  to  each 
other.  In  Figure  2.16  and  Figure  2.17,  the  shape  and  the  curvatures  at  the  end  of  the 
development  of  the  interface  are  plotted  while  the  former  shows  the  situation  before 
fairing  is  applied,  the  latter  illustrates  the  effect  after  the  fairing.  Furthermore,  with  a 
periodic  curve,  the  curvature  should  be  periodic  as  well.  As  shown  in  Figure  2.13(b), 
initially,  the  computed  curvature  is  periodic;  however,  without  fairing,  Figure  2.16(b) 
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Figure  2.12:   (a)  Shapes  of  expanding  circle  at  equal  intervals  of  time,    (b)  Corre- 
sponding curvatures  along  the  curve  at  the  same  time  instants  as  in  (a). 
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shows  that  the  computed  curvature  for  a  moved  interface  is  no  longer  periodic.  After 
fairing,  as  shown  in  Figure  2.17(b),  the  periodicity  is  restored  again. 

The  shape  is  well  preserved  and  the  oscillation  in  curvature  is  removed.  The 
majority  of  the  curvature  plot  after  fairing  is  the  same  as  that  before  fairing.  The 
number  of  marker  points  increases  as  a  result  of  reorganizing  process  to  follow  the 
increasing  arc  length  of  the  bulging  out  interface. 

The  B-spline  curve  fitting  in  conjunction  with  the  fairing  operation  is  a  robust 
interface  representing  and  tracking  method  which  can  yield  smooth  interface  and 
smooth  curvatures  while  preserve  the  interface  shape.  The  programming  task,  how- 
ever, is  quite  challenging.  For  closed  and  open  curves,  different  program  modules  and 
algorithms  have  to  be  implemented  for  B-spline  fitting  mainly  because  the  end  con- 
ditions are  different.  For  closed  curves,  the  end  condition  is  always  periodic  (cyclic). 
In  the  cases  of  open  curves,  there  might  be  other  end  conditions. 

2.3     Cartesian  Grid  Method  for  Sharp  Interfaces 

Once  the  interface  is  defined,  one  needs  to  identify  in  which  phase  each  compu- 
tational cell  lies  so  that  correct  transport  properties  can  be  assigned.  Furthermore, 
for  any  two  neighboring  cells  between  which  a  interface  passes  through,  there  may  be 
a  discontinuity  in  transport  properties  such  as  density  and  viscosity.  Those  boundary 
cells  are  reshaped  to  maintain  flux  conservation  around  the  interface.  A  major  goal 
of  the  present  approach  is  to  adopt  a  finite- volume  formulation  for  all  computational 
cells  so  that  mass,  momentum,  and  energy  conservation  is  honored  in  all  resolvable 
scale,  including  the  computational  cells  intersecting  with  phase  boundaries.  Once  a 
cell  intersects  with  an  interface,  it  is  split  into  two  parts;  with  the  partial  cell  con- 
taining the  center  of  the  original  Cartesian  cell  maintaining  the  initial  cell  index,  and 
the  other  merged  into  a  neighboring  cell  belong  to  the  same  phase.  As  illustrated 
in  Fig.  2.19,  this  procedure  results  in  irregular,  trapezoidal  shaped  cells  around  the 
interfaces.  Away  from  the  interfaces,  most  cells  are  still  structured  Cartesian  cells. 
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Figure  2.13:  (a)  The  initial  shape  of  two  fingers,  (b)  The  corresponding  curvature  of 
two  fingers  using  B-spline  fitting. 
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Figure  2.14:  (a)  The  shape  of  two  fingers  at  the  early  stage  of  the  development  when 
fairing  is  not  applied,  (b)  The  curvature  corresponding  to  the  shape  denoted  by  solid 
line  in  (a). 
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Figure  2.16:  (a)  The  shape  evolution  of  two  fingers  until  the  later  stage  of  the  de- 
velopment when  fairing  is  not  applied,  (b)  The  curvature  corresponding  to  the  final 
shape  denoted  by  solid  line  in  (a). 
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Figure  2.18:  (a)  The  corrections  of  the  interface  location  by  fairing  algorithm  at 
time  instants  corresponding  to  shapes  shown  in  2.17(a).  (b)  The  corrections  of  the 
interface  location  by  fairing  algorithm  at  the  same  markers  as  in  (a). 
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Figure  2.19:  Illustration  of  the  resulting  situation  when  cut-cell  approach  is  applied, 
i.e.  fragments  of  cells  which  are  cut  by  the  interface  are  absorbed  into  neighboring 
cells.  The  newly  formed  cells  are  shown  in  dashed  lines  on  both  sides  of  the  interface 


In  such  an  algorithm,  while  a  nominal  structured  grid  index  system  is  maintained, 
the  flux  calculation  across  the  surface  is  conducted  based  on  interpolation  of  varying 
degrees  of  polynomials.  The  resulting  flow  solver  needs  to  account  for  both  Cartesian 
and  trapezoid  cells,  using  a  finite  volume,  fractional  step  method. 

After  the  domain  is  partitioned  into  several  sub-domains  representing  different 
phases.  The  so-called  Cartesian  grid  method  is  developed  to  deal  with  the  discretiza- 
tion and  solution  of  governing  equations  in  such  a  domain  left  with  mostly  structured 
rectangular  cells  and  some  unstructured  trapezoidal  shape  cells  in  the  vicinity  of  the 
interface.  Essentially,  what  the  Cartesian  grid  method  does  is  to  solve  the  governing 
equations  on  a  domain  with  fixed,  complex  geometry  interfaces  without  resorting  to 
the  body-fitted  grid.  Although  the  method  is  developed  in  this  research  as  one  of 
the  major  components  in  the  entire  moving  interface  simulation  package,  the  method 
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alone  is  an  good  alternative  to  the  body-fitted  grid  approach  to  solve  the  flow  over 
solid  bodies  with  complex  geometry,  as  demonstrated  in  Ye  et  al.  (1999). 

This  globally  second-order  accurate  Cartesian  grid  flow  solver  for  fixed,  immersed 
interfaces,  (see  Ye  et  al.,  1999),  is  summarized  in  the  following: 

1.  The  interface  is  represented  numerically  using  C2  cubic  B-spline  interpolatory 
curve; 

2.  A  fractional  step  procedure  (see  Chorin,  1968;  Kim  &  Moin,  1985)  is  employed 
to  solve  the  coupled  governing  equations,  which  results  in  a  efficient  solution  of 
time  accurate,  unsteady  flows; 

3.  Colocated  variable  arrangement  (see  Ferziger  &  Peric,  1996)  is  adopted,  which 
leads  to  simplicity  in  coping  with  cut  cells  near  the  interface, 

4.  Second-order  accurate  central  difference  scheme  is  used  for  spatial  discretiza- 
tion, Crank-Nicolson  scheme  for  temporal  discretization; 

5.  A  compact  interpolation  scheme  near  the  immersed  interfaces  is  devised  that 
allows  us  to  retain  second-order  accuracy  and  conservation  property  of  the 
solver. 

6.  A  preconditioned  conjugate  gradient  iteration  method  is  used  for  solving  the 
Poisson  equation  for  pressure,  which  takes  advantage  of  the  underlying  struc- 
tured nature  of  the  mesh  and  also  substantially  accelerates  the  convergence  of 
the  Poisson  equation  for  pressure. 

2.3.1     Fractional  Step  Method  for  Unsteady  Transient  Solution 

To  use  finite  volume  discretization,  the  dimensionless  governing  equations  listed 
in  the  section  1.1.7  on  page  9  have  to  be  written  in  the  integral  forms  by  integrating 
each  term  over  the  control  volume  cv  and  wherever  possible,  using  Gauss  theorem 
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to  convert  volume  integrals  into  surface  integrals.  To  solve  the  momentum  equations 
(1.18)  and  (1.21)  under  the  constraints  (1.17)  and  (1.20)  in  the  liquid  and  vapor 
phases  respectively,  a  so  called  fractional  step  procedure  is  applied.  After  the  mo- 
mentum equations  are  solved,  the  energy  equations  (1.19)  and  (1.22)  can  then  be 
solved  in  the  updated  flow  field. 

All  governing  equations  are  discretized  on  a  Cartesian  mesh  using  a  colocated 
arrangement  (see  Ferziger  &  Peric,  1996)  of  the  primitive  variables  which  are  located 
at  the  cell  center.  For  convenience  of  describing  the  numerical  procedure,  only  the 
integral  form  of  dimensionless  governing  equations  in  the  liquid  phase  are  considered 
throughout  this  chapter  as  examples,  the  vapor  phase  transport  equations  are  solved 
in  the  exactly  same  procedure  except  additional  dimensionless  fluid  property  ratios 
are  present  on  those  cells  belong  to  the  vapor  phase. 

The  integral  form  of  Navier-Stokes  equations  in  the  liquid  phase  are: 

Mass  conservation: 


/ 


u-ndS  =  0  (2.10) 

cs 

Momentum  conservation: 


/  —  dV  +  I  u  (ti  •  n)  dS  =  -      pndS  +  —  /  Vu  •  ndS  (2.11 


In  the  above  equations  cv  and  cs  denote  the  control  volume  and  control  surface 
respectively  and  n  is  a  unit  vector  normal  to  the  control  surface.  The  above  equations 
are  to  be  solved  with  u(x,t)  =  v(x,t)  on  the  boundary  of  the  flow  domain  where  v 
is  the  prescribed  boundary  velocity. 

A  second  order  accurate,  two  step  fractional  step  method  (see  Chorin,  1968;  Kim 
&  Moin,  1985;  Zang  et  al.,  1994)  is  used  for  advancing  the  solution  in  time.  In  this 
time  stepping  scheme,  the  solution  is  advanced  from  time  step  'n'  to  'n+F  through 
an  intermediate  advection  diffusion  step  where  the  momentum  equations  without  the 
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pressure  gradient  terms  are  first  advanced  in  time.  The  semi-discrete  form  of  this 
advection  diffusion  equation  for  each  cell  shown  in  Figure  2.20  can  be  written  as: 

fU*~UndV    =    -)-f[?>un{Un-n)-un-l(Un-l-n)}dS 

cv  cs 

+r4-  /  (Vu*  +  Vun)  n  dS  (2.12) 

2Re  J 

cs 

where  u*  is  the  intermediate  cell-center  velocity.  A  second-order  Adams-Bashforth 
scheme  is  employed  for  the  convective  terms  and  the  diffusion  terms  are  discretized 
using  the  implicit  Crank-Nicolson  scheme.  This  eliminates  the  viscous  stability  con- 
straint which  can  be  quite  restrictive  in  the  simulation  of  viscous  flows. 

The  velocity  boundary  condition  imposed  at  this  intermediate  step  corresponds 
to  that  at  the  end  of  the  full  time  step,  i.e.,  u*  =  vn+l. 

At  this  stage,  in  addition  to  the  cell  center  velocities  which  are  denoted  by  u, 
we  also  introduce  face  center  velocities  U.  In  a  manner  similar  to  a  fully  staggered 
arrangement,  only  the  component  normal  to  the  cell  face  is  computed  and  stored 
(see  Figure  2.20).  The  face  center  velocity  is  used  for  computing  the  surface  flux 
from  each  cell  in  our  finite  volume  discretization  scheme.  The  reason  of  separately 
computing  the  face  center  velocities  will  be  addressed  later  in  this  section. 

Following  the  advection  diffusion  step,  the  intermediate  face  center  velocity  U* 
is  computed  by  interpolating  from  the  cell  center  intermediate  velocities. 

The  advection  diffusion  step  is  followed  by  the  pressure  correction  step 

un+l  -  u* 

£ -Vpn+1  (2.13) 

where  we  require  that  the  final  velocity  field  satisfy  the  integral  mass  conservation 
equation  given  by 


[  (Un+l-n)dS  =  0  (2.14) 
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Figure  2.20:  Schematic  of  underlying  Cartesian  mesh  and  arrangement  of  cell-center 
and  face-center  velocities 


This  results  in  the  following  equation  for  pressure 

j  (Vp)-ndS=±tJ  {U*-n)dS 


(2.15) 


which  is  the  integral  version  of  the  Poisson  equation  for  pressure.  Note  that  the 
pressure  correction  step  is  represented  by  the  inviscid  equation  (2.13)  and  is  well 
posed  only  if  the  velocity  component  normal  to  the  boundary  is  specified.  Therefore 
the  velocity  boundary  condition  consistent  with  Eq.  (2.13)  is  un+1-n  =  vn+l-n  where 
n  is  the  unit  normal  to  the  boundary  of  the  flow  domain.  It  can  be  easily  shown 
that  this  implies  that  (Apn+1)  n  =  0  be  used  as  the  boundary  condition  for  equation 
(2.15).  Once  the  pressure  is  obtained  by  solving  this  equation,  both  the  cell  center 
and  face  center  velocities  are  updated  separately  as  follows: 


un+1  =  u*  -  At  (Vpn+1)cc     ;     Un+l  =  U*  -  At  {Vpn+1)fc 


(2.16) 


63 

It  is  a  well  known  fact  that  in  a  colocated  mesh,  when  the  compact  stencil  is  used 
for  the  Poisson  equation  for  pressure,  the  pressure  field  usually  does  not  exhibit  non- 
physical  oscillation,  however,  the  final  velocity  field  does  not  satisfy  the  divergence 
free  constraint  exactly.  To  overcome  this  problem,  two  types  of  approaches  could  be 
used,  i.e.,  the  momentum  interpolation  (see  Rhie  &  Chow,  1983)  scheme  in  which  a 
pressure  gradient  term  is  introduced  into  the  interpolation  process  for  obtaining  cell 
face  velocities.  As  it  is  shown  by  Ferziger  &  Peric  (1996),  this  technique  produces  a 
third-order  dissipation  in  estimating  the  cell  face  velocities.  The  momentum  interpo- 
lation technique  is  not  straightforward  to  implement  in  cut  cells  encountered  in  the 
mixed  Eulerian-Lagrangian  algorithm  for  flow  with  immersed  interfaces. 

The  other  approach  for  dealing  with  this  problem  was  first  used  by  Zang  et  al. 
(1994).  In  their  method,  both  cell  center  and  face  center  velocities  are  stored  while 
momentum  equations  and  Poisson  equation  for  pressure  are  solved  at  the  cell  centers 
only.  After  Poisson  equation  for  pressure  is  solved,  both  cell  center  and  face  center 
velocities  are  updated.  The  pressure  gradient  at  the  cell  face  is  computed  in  a  compact 
stencil  manner.  They  have  used  this  procedure  in  conjunction  with  a  curvilinear  mesh 
solver  to  simulate  the  turbulent  flows  with  LES.  In  our  algorithm,  we  implemented 
this  method  and  test  results  show  that  this  approach  is  well  suited  for  high  Reynolds 
number,  unsteady  flows. 

2.3.2     Discretization  on  Regular  and  Cut  Cells 

To  discretize  governing  equations  in  the  configuration  as  shown  in  Figure  2.19, 
we  have  to  deal  with  two  types  of  situations:  discretization  on  cells  which  are  away 
from  the  interface  thus  not  cut  by  the  interface  and  on  those  cells  which  are  cut  by 
the  interface  and  reshaped  to  become  trapezoidal  cells. 

The  numerical  procedure  for  solving  Navier-Stokes  equations  just  outlined  can 
be  easily  implemented  on  a  regular  Cartesian  grid  cell.  In  this  section,  detailed 
finite  volume  integration  of  each  term  is  described.  For  demonstration  purpose,  basic 
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Figure  2.21:  Schematic  of  cells  around  the  interface,  (a)  boundary  cells  with  interface 
located  south  of  cell  center,  (b)  boundary  cells  with  interface  located  west  of  cell 
center,  (c)  typical  reshaped  quadrilateral  cells  corresponding  to  case  (a),  (d)  typical 
reshaped  quadrilateral  cells  corresponding  to  case  (b). 


procedure  for  evaluating  volume  and  surface  integrals  on  a  regular  cell  is  explained 
first,  followed  by  the  special  treatment  on  those  cells  cut  by  the  immersed  interfaces 
and  reshaped,  as  shown  in  Figure  2.19.   The  emphasis  is  placed  on  the  latter. 

In  general,  the  so  called  cut  cell  approach  is  used  to  treat  those  reshaped  cells: 
the  small  fragments  of  cut  cells  are  absorbed  into  neighboring  cells  in  the  same  phase 
to  form  mainly  trapezoidal  cells  as  shown  in  Figure  2.21,  after  that,  the  governing 
equations  are  integrated  on  those  reshaped  cells. 

It  is  straightforward  to  obtain  second  order  accurate  spatial  discretization  on 
regular  cells.  The  key  for  retaining  global  second-order  accuracy  of  the  solver  is  to 
evaluate  mass,  convective  and  diffusive  fluxes  and  pressure  gradients  to  second  order 
accuracy  on  the  cell  face  of  those  trapezoidal  cells  resulted  from  reshaping  process. 
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As  shown  in  equation  (2.12),  a  finite  volume  discretization  of  Navier-Stokes 
equations  requires  the  estimation  of  surface  integrals  on  the  faces  of  each  cell.  The 
integrand  (denoted  here  by  /  can  either  involve  the  value  or  the  normal  derivative  of 
a  variable.  An  example  of  the  former  is  the  convective  flux  denoted  by  (pcfrv-n)  and 
of  the  latter,  the  diffusive  flux  given  by  (rV</>-n)  where  <j>  is  a  generic  scalar  variable 
and  T  is  diffusive  coefficient.  In  addition  to  this,  the  Poisson  equation  for  pressure 
also  requires  evaluation  of  the  normal  pressure  gradient.  In  order  to  estimate  these 
surface  integrals  with  second  order  accuracy,  the  midpoint  rule  can  be  used  and  this 
requires  an  accurate  evaluation  of  the  integrand  at  the  center  of  the  face.  For  regular 
cells  which  are  away  from  the  interface  the  integrand  can  be  evaluated  at  the  face 
center  to  second  order  accuracy  in  a  straightforward  manner  by  assuming  a  linear 
profile  between  nodes  on  either  side  of  the  face. 

However,  this  is  not  the  case  for  the  trapezoidal  cells.  Unlike  regular  Cartesian 
cells,  for  those  trapezoidal  cells  near  the  interface,  simple  linear  approximation  be- 
tween neighboring  cell  centers  can  not  be  made  to  calculate  the  fluxes  and  pressure 
gradients  on  the  face  centers  to  second  order  accuracy.  This  is  because  the  centers 
of  some  of  the  faces  of  such  a  cell  (marked  by  a  shaded  arrow  in  Figure  2.21c  and 
d)  may  not  lie  in  a  location  which  puts  it  in  the  middle  of  neighboring  cell  centers 
where  a  linear  approximation  would  produce  second  order  accurate  estimate  of  the 
gradients.  Furthermore,  some  of  the  neighboring  cell  centers  do  not  even  lie  on  the 
same  side  of  the  immersed  interface  and  therefore  cannot  be  used  in  the  differencing 
procedure.  Thus,  not  only  do  we  need  a  procedure  for  computing  these  face  center 
quantities  accurately,  but  we  also  require  that  the  procedure  adopted  be  capable  of 
systematically  handling  reshaped  boundary  cells  with  a  wide  range  of  shapes. 

Our  solution  has  been  to  use  a  compact  two  dimensional  polynomial  interpo- 
lating function  which  allows  us  to  obtain  a  second  order  accurate  approximation  of 
the  fluxes  and  gradients  on  the  faces  of  the  trapezoidal  boundary  cells  from  available 
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Figure  2.22:  Schematic  of  interpolation  for  cell  face  values  and  derivatives  at  trape- 
zoidal cells  (a)  various  fluxes  required  for  trapezoidal  cell  (b)  trapezoidal  region  and 
stencil  used  in  computing  fsw. 


neighboring  cell  center  values.  The  current  interpolation  scheme  coupled  with  the 
finite  volume  formulation  guarantees  that  the  accuracy  and  conservation  property  of 
the  underlying  algorithm  are  retained  even  in  the  presence  of  curved  immersed  inter- 
faces. In  the  following,  we  describe  the  interpolation  function  for  a  typical  trapezoidal 
boundary  cell. 

Consider  the  trapezoidal  cell  ABCVS  in  Figure  2.22a.  The  face  ABC  of  the 
trapezoidal  cell  is  composed  of  two  pieces:  AB  coming  from  cell  P  and  BC  coming 
from  cell  S.  The  integral  on  this  face  can  be  decomposed  as 

f  fdy  =  j  fdy  +  J  fdy  (2.17) 

AC  AB  BC 

A  second  order  approximation  to  this  integral  can  then  be  obtained  as 


/ 

AC 


fdy&fw  (yA  -  Vb)  +  fsw  {Vb  -  Vc) 


(2.18) 


where  /wand  fsw  are  computed  at  the  centers  of  segments  AB  and  BC  respectively.  On 
the  other  hand,  if  the  face  is  cut  by  the  immersed  interface  such  that  it  is  smaller  than 
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a  nominal  cell  face,  as  in  the  case  of  face  VS  then  the  integral  can  be  approximated 
as 


/ 


fdy*fe(y£-yv)  (2.19) 

vs 
where  fe  is  the  flux  computed  at  the  center  of  the  segment  V£.  For  non-boundary 

cells,  these  face  center  values  can  be  evaluated  to  second  order  accuracy  quite  easily 

by  a  linear  approximation  and  we  would  therefore  like  to  evaluate  /„,,  fsw,  fe  to  within 

second  order  accuracy  also.  Approximation  of  fw  to  second  order  accuracy  is  quite 

straightforward  and  is  done  in  the  same  way  as  for  the  face  of  a  non-boundary  cell. 

For  instance,  if  /„,  requires  the  value  of  </>,  this  can  be  evaluated  to  second  order 

accuracy  as 

<f>w  =  <f>w^w  +  4>p  (1  -  ^w)  (2.20) 

where  the  linear  interpolation  factor  Aw  is  defined  as 

Xw  =  Xp~Xyf  (2.21) 

Xp  —  X\y 

Alternatively,  if  fw  requires  the  normal  gradient  of  </»  as  it  would  for  the  diffusion  or 
pressure  gradient  terms,  this  can  be  approximated  by  a  central  difference  scheme  as 
follows: 

(*£\     =^~^  (2.22) 

\dxjw      xp-xw 

This  approximation  is  second  order  accurate  when  the  cell  face  is  midway  between 
P  and  W,  i.e.,  when  the  mesh  is  uniform.  However,  expressions  such  as  Eq.  (2.20) 
or  Eq.  (2.22)  cannot  be  used  in  evaluation  of  fsw  or  fe  to  second  order  accuracy 
since  in  many  instances  some  of  the  neighboring  nodes  can  be  inside  the  immersed 
interface.  For  instance,  for  the  situation  shown  in  Figure  2.22a,  the  south  node  is 
inside  the  immersed  interface  and  cannot  be  used  in  the  evaluation  of  fsw.  Even 
if  neighboring  nodes  are  available,  as  they  are  for  the  east  face,  it  is  not  clear  how 
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a  second  order  accurate  scheme  can  be  constructed  since  fe  is  not  located  on  the 
line  joining  the  neighboring  cell  centers  and  consequently,  expressions  such  as  Eq. 
(2.20)  or  Eq.  (2.22)  cannot  approximate  this  flux  to  second  order  accuracy.  Thus,  a 
different  approach  is  needed  here  for  evaluating  these  fluxes. 

Our  approach  is  to  express  the  flow  variables  in  terms  of  a  two  dimensional 
polynomial  interpolating  function  in  an  appropriate  region  and  evaluate  the  fluxes 
such  as  fsw  or  fe  based  on  this  interpolating  function.  For  instance,  in  order  to 
approximate  fsw,  we  express  4>  m  the  shaded  trapezoidal  region  shown  in  Figure 
2.22b  in  terms  of  a  function  that  is  linear  in  x  and  and  quadratic  in  y 

<f>  =  cixy2  +  c2y2  +  c3xy  +  c4y  +  c5x  +  c6  (2.23) 

where  ex  to  ce  are  six  unknown  coefficients.  If  fsw  involves  the  normal  derivative  of 
(j),  this  can  be  obtained  by  differentiating  the  interpolating  function,  i.e. 

^  =  clV2  +  c3y  +  c5  (2.24) 

The  rationale  for  choosing  Eq.  (2.23)  as  the  interpolating  function  for  evaluating 
faw  is  as  follows:  the  objective  here  is  to  evaluate  (d<f)/dx)  at  the  center  of  BC  to 
within  at  least  second  order  accuracy.  Furthermore,  we  would  like  to  do  this  with 
the  most  compact  interpolating  function  so  as  to  minimize  the  size  of  the  stencil 
required  for  the  boundary  cell.  Clearly,  a  biquadratic  interpolating  function  in  the 
trapezoid  shown  in  Figure  2.22b  would  lead  to  a  second  order  accurate  evaluation  of 
the  derivative  anywhere  inside  the  trapezoid.  However,  a  biquadratic  function  has 
nine  unknown  coefficients  and  therefore  requires  a  large  nine  point  stencil.  It  turns 
out  however  that  for  the  trapezoid  shown  in  Figure  2.22b,  second  order  accurate 
evaluation  of  the  derivative  on  the  cell  face  can  be  achieved  by  using  an  interpolating 
function  that  is  quadratic  in  y  but  only  linear  in  x.  This  is  because  BC  is  midway 
between  the  two  parallel  sides  of  the  trapezoidal  and  in  a  manner  analogous  to  central 
differencing,  linear  interpolation  in  the  ^-direction  leads  to  a  second  order  accurate 
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evaluation  of  derivative  at  this  location.  On  the  other  hand,  this  situation  does  not 
exist  in  the  y-direction  for  the  cell  shown  in  Figure  2.22b  and  therefore  a  quadratic 
interpolation  is  necessary  in  this  direction  in  order  to  obtain  a  second-order  accurate 
approximation  to  (d<fi/dx)  at  the  center  of  BC.  In  Appendix  1,  we  have  demonstrated 
numerically  that  the  linear  quadratic  interpolating  function  shown  in  Eq.  (2.23)  does 
indeed  result  in  second  order  accurate  evaluation  of  values  and  derivatives  on  a  line 
which  is  located  midway  between  the  two  parallel  sides  of  a  trapezoid. 

It  can  be  seen  in  Figure  2.22b  that  the  sides  of  the  trapezoid  in  which  the 
interpolation  is  performed  pass  through  four  nodal  points  and  two  boundary  points. 
Thus,  the  six  unknown  coefficients  in  2.23  can  be  expressed  in  terms  of  the  values  of 
<f>  at  these  six  locations.  To  solve  for  cn,  we  obtain  the  following  system  of  equation 
by  expressing  the  </>  at  the  six  location  in  terms  of  the  linear  quadratic  interpolating 
function: 


^2  y2    y\    x2y2    y-i     x2     1 


c2 
;  ce  J 


(2.25) 


XsVl      Vl      x62/6      2/6       x6        1 

The  coefficients  can  now  be  expressed  in  terms  of  values  of  <j)  at  the  six  points  by 
inverting  Eq.  (2.25),  i.e. 


cn  =  ^bnjfij,     n  =  1,6 
j=i 


(2.26) 


where  bnj  are  the  elements  of  the  inverse  of  the  matrix  in  Eq.  (2.25). 

After  cn  is  obtained,  the  value  of  <f>  at  the  center  of  BC  is  expressed  in  the  form 
of 


<f>sw  =  cixswylw  +  c2y2sw  +  c3xswysw  +  c4ysw  +  cbxsw  +  c6 
and  using  Eq.  (2.26)  this  can  be  rewritten  as 

6 
<f>sw  =  22  Oij<f)j 


(2.27) 


(2.28) 
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where 

Ctj  -  hjXawVtm  +  b2jV2sw  +  hjXswVsw  +  hjVsw  +  hj^sw  +  hj  (2.29) 

The  value  of  (d<p/dx)  at  the  center  of  BC  is  expressed  as 

[a-)      =  c^vL  +  c*Vsw  +  c5  (2.30) 


sw 


6 


and  using  Eq.  (2.26)  this  can  be  written  as 

B)„-'S**  (2'31) 

where 

Pj  =  Kvlw  +  hjVsw  +  hj  (2.32) 

Note  that  a  and  j3  are  coefficients  that  depend  only  on  the  mesh,  the  location  and 
orientation  of  the  immersed  interface.  Therefore  these  can  be  computed  once  at  the 
beginning  of  the  solution  procedure.  Subsequently,  relationships  such  as  (2.28)  and 
(2.31)  can  be  used  in  the  spatial  discretization  of  the  governing  equations  (2.12)- 
(2.16). 

A  similar  interpolation  procedure  is  also  used  for  approximating  fe.  For  this, 
a  linear  quadratic  interpolating  function  is  used  in  the  trapezoidal  region  shown  in 
Figure  2.23  and  a  relationship  similar  to  Eq.  (2.28)  and  (2.31)  is  developed  for  ap- 
proximating fe.  The  six  points  contained  in  this  stencil  are  also  shown  in  Figure 
2.23.  It  should  be  pointed  out  that  the  north  face  of  the  particular  cell  being  con- 
sidered here  does  not  need  special  treatment  since  face  center  values  and  derivatives 
can  be  computed  to  second  order  accuracy  using  a  linear  approximation.  However, 
in  general  there  are  also  boundary  cells  which  have  their  north  or  south  faces  cut  by 
the  immersed  interface  (as  shown  in  Figure  2.21b).  For  these  boundary  cells  too,  the 
same  approach  is  used  to  evaluate  the  fluxes  on  these  cut  faces.  The  only  difference 
here  is  that  the  interpolating  function  is  linear  in  y  and  quadratic  in  x. 
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Figure  2.23:  Trapezoidal  region  and  stencil  used  in  computing  fe 


Now  we  turn  to  the  calculation  of  the  flux  on  cell  face  CV  which  lies  on  the 
immersed  interface  as  shown  in  Figure  2.22a.  The  integrated  flux  on  this  face  can 
again  be  evaluated  to  second  order  accuracy  using  the  midpoint  rule  and  as  before 
we  would  like  to  evaluate  the  integrand  at  the  center  of  face  CV  (denoted  here  by 
fint)  to  second  order  accuracy.  In  general  both  convective  and  diffusive  fluxes  are 
needed  on  this  face  and  this  requires  approximation  of  variable  values  as  well  normal 
derivatives  at  the  center  of  CV.  The  value  is  usually  available  from  a  Dirichlet  type 
boundary  condition  and  hence  no  interpolation  is  required  for  this.  Here  we  describe 
the  approximation  procedure  for  the  normal  derivative.  The  normal  derivative  on 
face  CV  can  be  decomposed  as 


d(j)      d(j)  d<f) 

dn      dx  dy    y 


(2.33) 


where  nx  and  ny  are  the  two  components  of  the  unit  vector  normal  to  face  CV. 
Since  we  know  the  shape  of  the  immersed  interface,  nx  and  ny  are  known.  Therefore 
computation  of  the  normal  flux  requires  estimation  of  d<j)/dx  and  d<j)/dy  at  the  center 
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Figure  2.24:  Stencil  for  calculating  interface  flux,   (a)  stencil  for  calculating  d<fi/dy 
(b)  stencil  for  calculating  dtfr/dx 


of  the  line  segment  CV.  For  the  cell  being  considered  here,  dcj)/dy  is  computed  to 
second  order  accuracy  with  relative  ease  by  expressing  the  4>  variation  along  the 
vertical  line  in  terms  of  a  quadratic  in  y  as  follows 


cf)  =  axy  +  a2y  +  a3 


(2.34) 


The  coefficients  in  the  quadratic  can  be  expressed  in  terms  of  the  values  of  t\>  at  the 
three  points  indicated  in  Figure  2.24a.  Subsequently,  the  normal  derivative  at  the 
center  of  face  CV  is  evaluated  as: 


\dy)     =2°1^  +  a2  =  X]rJ'^ 


(2.35) 


where  again  rj  are  coefficients  which  depend  solely  on  the  geometry  of  the  boundary 
cell. 

Unlike  the  calculation  of  d<p/dy  for  this  cell,  the  calculation  of  d<j)/dx  is  not 
straightforward.  However,  an  approach  consistent  with  the  computation  of  fsw  and 
/e  can  be  used  to  estimate  the  value  of  this  derivative  to  desired  accuracy.  Consider 
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the  trapezoid  shown  in  Figure  2.24b.  Again,  because  the  x-coordinate  of  the  center 
of  CV  is  located  midway  between  the  two  parallel  sides  of  this  trapezoid,  expressing 
the  variable  in  this  trapezoid  in  terms  of  an  interpolating  function  which  is  linear  in 
x  and  quadratic  in  y  allows  us  to  obtain  a  second  order  accurate  approximation  to 
(d(j)/dx)sw  at  the  center  of  the  line  segment  CV.  The  procedure  for  this  follows  along 
lines  similar  to  that  shown  earlier  for  (d<j)/dx)sw  and  we  get  the  following  expression 
for  the  x-derivative  on  the  interface: 

where  if  depend  on  the  location  and  orientation  of  the  immersed  interface  in  the 
neighborhood  of  the  cell  under  consideration.  Finally  we  obtain  an  expression  of  the 
form 

6EL-g* 

for  the  normal  gradient  where  Tj  can  be  obtained  from  Eq.  (2.33),  (2.35)  and  (2.36). 
Thus  we  obtain  a  nine-point  stencil  for  the  flux  on  the  interface  and  the  points  in  this 
stencil  are  shown  in  Figures  2.24a  and  b.  As  can  be  seen  from  these  figures,  of  these 
nine  points,  three  points  lie  on  the  immersed  interface  and  their  values  are  available 
from  the  prescribed  boundary  condition. 

It  should  be  pointed  out  that  although  most  cells  are  four-sided  trapezoidal  cells, 
some  five-sided  cells  and  three-sided  triangular  cells  may  also  be  encountered.  How- 
ever, the  discretization  of  the  governing  equations  for  these  cells  can  also  be  handled 
within  the  framework  of  the  current  interpolation  scheme.  With  all  of  these  features, 
the  current  solver  can  in  principle,  handle  arbitrarily  complex  geometries.  Further- 
more, multiple  immersed  interfaces  can  be  handled  as  easily  as  a  single  interface. 
This  is  in  contrast  to  body  fitted  grid  where  the  grid  topology  can  get  quite  compli- 
cated in  the  presence  of  multiple  interfaces.  Finally,  since  both  sides  of  the  immersed 
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interface  are  covered  by  the  grid,  we  have  the  capability  to  solve  a  different  set  of 
equations  on  each  side  of  the  immersed  interface.  For  instance,  the  multi-phase  flow 
with  separate  governing  equations  in  different  phases  across  the  interface  could  be 
simulated  in  this  manner. 

2.3.3     Imposition  of  Boundary  Conditions  on  the  Interface 

This  section  describes  the  imposition  of  boundary  conditions  on  curved,  immersed 
interfaces  on  the  underlying  Cartesian  grid 

In  the  discretization  process  on  cut  cells  as  discussed  in  the  last  section,  the 
interpolation  scheme  used  to  obtain  the  face  flux  such  as  fe,  fsw  and  fint  as  shown  in 
Figure  2.22  involves  the  values  of  generic  variable  (j)  at  point  (5),  (6)  on  the  interface 
as  illustrated  in  Figure  2.22,  2.23  and  2.24.  The  values  of  (j)  on  the  interface  are 
simply  known  if  Dirichlet  boundary  conditions  are  specified  on  the  interface.  In  the 
case  of  Neumann  boundary  conditions  on  the  interface,  the  explicit  values  of  <j>  on 
the  interface  need  to  be  extrapolated  from  values  in  the  interior  cells. 

The  extrapolation  process  is  explained  in  Figure  2.25.  For  implementing  Neu- 
mann boundary  conditions  on  the  interface,  it  is  necessary  to  obtain  the  expression 
for  gradients  |£  on  the  interface.  Taking  one  side  of  the  interface  as  example  as  shown 
in  Figure  2.25,  the  gradient  can  be  evaluated  by  using  a  normal  probe,  i.e.,  drawing  a 
normal  into  one  phase  from  the  interfacial  marker  point  location.  Since  second  order 
accuracy  is  desired,  two  nodes  on  the  normal  probe  are  needed  which  are  located  at 
distance  h  and  2/i  from  the  interface.  Values  at  each  node,  0nl|  <j)n2  are  obtained  by 
bilinear  interpolation  from  the  neighboring  grid  points  which  are  denoted  by  shaded 
circles.  Then  based  on  the  values  at  two  points  on  the  normal  probe  and  the  known 
gradient  values  on  the  interface,  an  0(h2)  estimate  of  the  normal  gradient  at  the 
marker  point  k  location  can  be  written  as: 

On       ~  2h  (2-38) 


75 


2nd  normal  probe  "*. 
node  <Pn2 

1st  normal  probc*^ 
node0„1 

O 

o 

Marker  points  k 

o 

o 

/ 

& 

r 

—  T. 

JX"*^ 

\f> 

_/ 

\ 

X  ^ 

Figure  2.25:  Illustration  of  implementing  Neumann  boundary  conditions  on  the  in- 
terface. 


Suppose  the  Neumann  boundary  condition  on  the  interface  gives  (f^) .  =  0,  the  value 
of  <f>  at  marker  point  k  can  be  extrapolated  based  on  this  expression  for  gradient  as 

4<^nl  -<j>n2-0-2h 


<f>a  = 


(2.39) 


At  each  marker  point,  the  value  4>  is  extrapolated  in  the  same  manner.  All  values  <j) 
at  sequential  marker  points  are  parameterized  in  terms  of  arc  length  s  as 


$(s)  =  a^s2  +  b^s  +  fy 


(2.40) 


In  the  cut  cell  procedure  described  earlier,  however,  the  values  of  <j>  at  the  inter- 
sections of  the  interface  with  grid  lines  are  needed,  e.g.,  05i6  in  the  Figure  2.25. 

Let  [xitJJk)  be  the  coordinate  of  this  interfacial  marker  point.  Previously  the 
functions  x(s)  and  y(s)  on  the  interface  for  each  marker  point  has  been  obtained. 
Thus,  the  arclength  value  s5>6  along  the  interface  where  the  curve  intersects  the  line 
x  —  Xi  is  obtained  by  solving 


(ax)kslfi  +  (bx)ks5fi  +  (cx)k  =  Xi 


(2.41) 
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Once  s5)6  is  obtained,  based  on  the  parameterization  of  <f>  in  terms  of  s,  the  value  of 
^5,6  is 

(f>5,6  =  Mk4,6  +  (MfcS5,6  +  (c^)fc  (2.42) 

These  values  can  then  be  incorporated  into  the  interpolation  scheme  to  calculate 
the  flux  on  the  cut  cell  faces  in  the  discretization  process. 

2.3.4     Inversion  of  Discrete  Operators 

This  section  discusses  our  choice  of  algebraic  equations  solver. 
The  finite  volume  discretization  of  Eq.  (2.12)  or  (2.15)  in  a  given  cell  P  can  be 
written  as 

M 

Y^Xptk  =  bP  (2.43) 

k=l 

where  x's  denote  the  coefficients  of  the  nodal  values  within  a  stencil  of  size  M  and 
bP  is  the  source  term  that  contains  the  explicit  terms  as  well  as  the  terms  involv- 
ing boundary  conditions.  The  term  on  the  left-hand  side  represents  a  discretized 
Helmholtz  operator  in  the  case  of  the  advection-diffusion  equation  and  a  Laplacian 
operator  in  the  case  of  the  Poisson  equation  for  pressure.  In  the  present  method,  M 
is  equal  to  five  for  non-boundary  cells  and  this  five-point  stencil  is  shown  in  Figure 
2.20.  For  the  trapezoidal  boundary  cells,  the  linear  quadratic  interpolation  scheme 
described  earlier  results  in  a  six-point  stencil  and  for  the  particular  boundary  cell 
considered  in  Figure  2.22,  the  resultant  six  point  stencil  is  shown  in  Figure  2.26. 

The  discretized  advection-diffusion  equation  and  Poisson  equation  for  pressure 
result  in  a  coupled  system  of  linear  algebraic  equations  which  requires  the  inversion  of 
a  large,  sparse,  banded  matrix.  The  structure  of  this  matrix  is  for  the  most  part  simi- 
lar to  that  obtained  on  a  Cartesian  mesh  without  any  immersed  interfaces.  However, 
the  presence  of  the  immersed  interface  alters  this  matrix  since  the  coefficients  in  Eq. 
(2.43)  are  different  for  the  trapezoidal  boundary  cells.  Furthermore,  the  rows  in  the 
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Figure  2.26:  Six  points  stencil  of  finite  volume  integration  on  a  typical  boundary  cell 


matrix  corresponding  to  the  trapezoidal  boundary  cells  also  have  additional  elements 
since  the  stencil  of  the  trapezoidal  boundary  cell  is  different  from  regular  cells.  The 
alternative  direction  line  successive-overrelaxation  (SOR)  method  (see  Press  et  a/., 
1992)  is  one  of  the  most  widely  used  iterative  methods  for  solving  equations  resulting 
from  discretization  on  structured  grids.  A  typical  implementation  of  this  technique, 
which  we  have  adopted  here,  involves  alternating  sweeps  along  the  two  families  of 
grid-lines.  We  find  that  even  in  the  presence  of  the  immersed  interfaces,  this  method 
is  extremely  effective  for  the  numerical  solution  of  the  discretized  advection-diffusion 
equation  and  the  residual  can  be  reduced  to  an  acceptable  level  within  a  few  itera- 
tions. 

The  discretized  Poisson  equation  for  pressure  however  exhibits  a  slower  conver- 
gence than  the  advection-diffusion  equation.  This  is  because  the  time-derivative  term 
in  the  advection-diffusion  equation  tends  to  improve  the  diagonal  dominance  of  the 
corresponding  discretized  operator.  In  fact  due  to  its  slow  convergence,  the  solution 
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of  the  discretized  Poisson  equation  for  pressure  is  usually  the  most  time  consum- 
ing part  of  a  fractional-step  algorithm.  In  the  presence  of  immersed  interfaces  this 
behavior  of  the  Poisson  equation  for  pressure  can  be  further  exacerbated  since  the 
stencil  for  the  trapezoidal  cells  contains  dependence  on  some  neighboring  cells  which 
are  not  included  in  the  line-SOR  sweeps.  For  example,  in  Figure  2.26  the  coupling  of 
cell  P  with  cells  (1)  and  (4)  in  the  stencil  is  not  included  directly  in  any  of  the  line 
sweeps.  Furthermore,  depending  on  the  aspect  ratio  of  the  trapezoidal  boundary  cell 
and  the  angle  at  which  the  immersed  interface  cuts  the  cell,  diagonal  dominance  in 
the  pressure  operator  can  be  severely  weakened.  In  the  various  simulations  that  have 
been  performed  using  the  current  method,  it  has  been  found  that  the  simple  line-SOR 
procedure  can  result  in  an  extremely  slow  convergence  for  the  pressure  equation. 

One  remedy  in  such  a  situation  is  to  resort  to  a  multigrid  method  (see  Brandt, 
1977;  Briggs,  1987).  However,  the  presence  of  the  immersed  interface  can  complicate 
the  implementation  of  a  multigrid  procedure  since  operations  such  as  prolongation 
and  restriction  are  difficult  to  perform  near  the  immersed  interface.  In  contrast, 
Krylov  subspace  methods  (see  Golub  &  Loan,  1989)  are  an  attractive  alternative 
since  these  are  designed  for  general  sparse  matrices  and  therefore  do  not  assume  any 
structure  in  the  matrix  operator.  Thus,  the  presence  of  the  immersed  interface  does 
not  pose  any  additional  complication  for  the  implementation  of  these  methods.  A 
particularly  suitable  method  in  this  class  is  the  bi-conjugate  gradient  stabilized  (Bi- 
CGSTAB)  method  (see  van  der  Vorst,  1992;  Barrett  et  a/.,  1993)  which  is  applicable 
to  non-symmetric  matrices  and  provides  relatively  uniform  convergence.  However,  as 
with  all  conjugate  gradient  methods,  the  convergence  rate  of  the  Bi-CGSTAB  proce- 
dure depends  critically  on  the  choice  of  the  preconditioner.  The  Jacobi  preconditioner 
which  has  a  trivial  construction  phase  is  used  routinely  in  conjunction  with  unstruc- 
tured grids.  However,  this  preconditioner  only  produces  marginal  improvement  in 
the  convergence  rate  of  conjugate  gradient  type  algorithms.   Other  preconditioners 
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such  as  those  based  on  incomplete  factorization  can  substantially  increase  the  con- 
vergence rate  but  these  usually  have  a  non-trivial  and  expensive  construction  phase 
(see  Barrett  et  al,  1993). 

The  structure  of  the  matrix  that  results  from  the  current  Cartesian  grid  ap- 
proach however  presents  us  with  another  alternative  choice  of  a  preconditioner.  As 
pointed  out  before,  the  presence  of  an  immersed  interface  only  alters  the  underlying 
matrix  operator  in  the  rows  corresponding  to  the  boundary  cells.  Although  this  alter- 
ation slows  the  convergence  rate  of  the  line-SOR  procedure,  the  convergence  rate  is 
still  significantly  faster  than  what  is  obtained  with  a  simple  point  Jacobi  method.  It 
follows  that  the  line-SOR  procedure  would  serve  as  a  better  preconditioner  than  the 
Jacobi  preconditioner.  A  further  advantage  of  using  the  line-SOR  as  a  preconditioner 
is  that  this  procedure  only  requires  the  solution  of  tridiagonal  systems  and  this  can  be 
accomplished  with  ease  using  the  Thomas  algorithm.  Thus,  in  the  solver  developed 
here,  the  line-SOR  procedure  was  used  as  a  preconditioner  in  the  Bi-CGSTAB  algo- 
rithm and  a  significant  improvement  in  the  convergence  rate  over  a  simple  line-SOR 
iterative  procedure  was  found. 

This  completes  the  discretization  and  solution  procedure  of  Navier-Stokes  equa- 
tions in  a  domain  with  immersed,  stationary  interfaces.  The  same  solution  procedure 
can  be  applied  to  solve  energy  equations  considering  it  is  an  advection-diffusion  equa- 
tion similar  to  momentum  equations. 

When  the  immersed  interfaces  are  no  longer  stationary  due  to  various  physical 
mechanisms  such  as  phase  change  or  hydrodynamic  forces  effect,  the  situation  be- 
comes much  more  complicated.  Many  issues  associated  with  solving  moving  interface 
problems  arise.  If  the  discontinuity  of  fluid  properties  become  appreciable  and  even 
severe,  the  overall  numerical  procedure  will  face  additional  challenges. 

In  the  following  sections,  key  issues  regarding  solving  transport  equations  for 
problems  with  moving  interfaces  are  addressed. 
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2.4     Moving  Interface  Algorithms 

In  the  previous  chapters,  the  numerical  procedure  used  to  solve  the  coupled, 
nonlinear  transportation  equations  in  the  domain  with  immersed  fixed  interfaces  is 
described.  The  method  is  fixed-grid  approach  in  which  the  grid  lines  don't  conform 
to  the  interface  with  complex  geometry. 

When  the  immersed  interfaces  in  the  flow  domain  translate  over  the  underlying 
fixed  Cartesian  grid  due  to  physical  mechanisms,  various  computational  issues  related 
to  the  moving  interface  can  arise. 

Overall,  the  numerical  procedure  involving  immersed  moving  interfaces  on  a  fixed 
Cartesian  grid  is  to  solve  a  fixed  interface  problem  at  any  time  instant  provided  that 
the  instantaneous  interface  location  is  known.  However,  since  the  interface  movement 
is  not  independent  of  the  flow  field,  this  coupling  between  the  interface  evolution  and 
the  flow  field  makes  the  moving  interface  problem  not  a  simple  addition  of  many 
individual  fixed  interface  problems  because  determining  the  interface  location  at  any 
time  step  is  also  part  of  the  solution  process.  For  fixed  interface  problems,  the  global 
iteration  process  is  only  carried  out  among  the  coupled  mass,  momentum  and  energy 
equations  to  solve  for  the  flow  field.  With  the  introduction  of  moving  interfaces,  the 
interface  location  is  coupled  with  the  flow  field,  another  outer  iteration  is  needed  to 
deal  with  this  higher  level  of  coupling.  The  major  difficulty  of  the  entire  simulation 
methodology  arises  in  this  iteration  process  to  solve  for  the  flow  field  and  to  adjust 
the  interface  shape  together. 

Physically,  the  interface  location  is  to  be  determined  based  on,  in  the  case  of 
liquid-vapor  two  phase  flow  considered  here,  the  thermal  effects  such  as  evaporation 
or  condensation  of  liquid  or  vapor  as  well  as  hydrodynamic  forces  such  as  pressure, 
viscous  stresses  in  the  two  phases.  The  interface  must  be  advected  to  a  position  where 
multiple  interfacial  conditions  such  as  mass,  momentum  and  energy  conservation 
conditions  (see  Equations  1.23,  1.24,  1.25),  and  mass  conservation  constraint  within 
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each  phase  are  simultaneously  satisfied.  In  order  to  achieve  this,  a  major  numerical 
technique  is  required  to  determine  the  interface  location  in  a  iterative  manner  within 
any  particular  time  step. 

2.4.1     Change  of  Phase 

As  it  is  mentioned  earlier  in  this  chapter,  when  a  interface  which  represents 
the  phase  front  of  the  liquid  and  vapor  state  moves  across  the  underlying  stationary 
Cartesian  grid,  some  grid  points  may  undergo  discontinuous  change  of  the  phase  due 
to  the  interface  motion  as  illustrated  in  the  Figure  2.27.  At  the  time  level  n,  the 
interface  is  located  such  that  the  grid  point  (i,j)  lies  in  the  vapor  phase  as  shown  in 
Figure  2.27a.  At  the  next  time  level  n+1,  the  interface  moves  to  the  updated  location 
shown  in  Figure  2.27b,  as  can  be  seen,  some  grid  points  which  previously  lie  in  the 
vapor  phase  now  lie  in  the  liquid  phase.  In  the  time  advancement  solution  procedure 
for  the  transportation  equations,  the  value  of  any  variable  4>  at  previous  time  level  n 
is  needed  at  current  time  level  n  +  1  by  the  temporal  discretization  scheme  adopted. 
However,  for  grid  point  that  have  just  changed  the  phase  since  last  time  level  n  such 
as  (i,  j)  in  the  Figure  2.27,  the  current  computation  at  time  level  n  +  1  is  performed 
in  the  liquid  phase  while  previous  history  for  this  point  is  in  the  vapor  phase  which 
couldn't  be  used,  i.e.,  there  is  no  physically  meaningful  information  regarding  the 
previous  value  of  <fin  at  (i,j).  This  is  a  purely  computational  issue  encountered  in  the 
numerical  approach  that  treats  the  interface  as  sharp  discontinuity  in  terms  of  the 
fluid  properties.  In  the  purely  Eulerian  methods,  the  interface  is  not  explicitly  tracked 
but  deduced  from  certain  field  variables  which  are  smoothed  over  the  fixed  grid,  thus 
the  interface  is  spread  over  a  band  of  grid  points  rather  a  sharp  discontinuity.  In  that 
case,  the  problem  of  change  of  the  phase  is  not  a  significant  issue  if  it  is  a  issue  at 
all.  However,  if  the  sharp  discontinuity  of  the  interface  is  desired  to  be  maintained 
like  in  our  approach,  those  grid  points  have  to  be  dealt  with. 
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A  solution  to  this  problem  is  to  estimate  the  value  of  any  variable  0  at  the  grid 
point  («,  j)  in  the  Figure  2.27  by  interpolation.  At  time  level  n  +  1,  the  grid  point 
(t,  j)  falls  into  liquid  phase,  but  this  point  was  in  the  vapor  phase  at  last  time  level  n. 
Therefore  the  value  of  tit*  doesn't  exist  in  terms  of  the  liquid  phase  and  is  obtained 
through  interpolation  between  the  surrounding  grid  points  in  the  liquid  phase  and 
points  on  the  interface. 

The  procedure  of  the  interpolation  is  as  following.  First  a  normal  to  the  interface, 
which  passes  through  the  grid  point  (i,j)  and  has  a  length  of  the  grid  spacing  h,  is 
drawn.  The  end  of  the  normal  is  called  /3,  the  intersection  of  the  normal  with  the 
interface  is  denoted  as  a.  The  value  of  0"  •  is  estimated  by  interpolation  through 
values  at  a  and  f3  as 


flj  =  a0a  + 

where  a,  b  are  interpolation  coefficients  based  on  length  ratios  between  a,  (3  and  (i,  j). 
4>a  is  the  known  boundary  value  on  the  interface,  <f>p  could  be  obtained  by  constructing 
a  bilinear  interpolation  within  point  1,  2, 3, 4  illustrated  in  the  Figure  2.27: 

£l#J  +  602  +  603  +  604  =  08 

where  the  fi-4  are  geometric  coefficients  in  the  bilinear  interpolating  function.  Com- 
bining these  two  expressions  gives  the  estimated  value  at  the  grid  point  (i,  j)  as 

. n       _  Q0Q  +  b  (602  +  603  +  £404) 

^  "  1-66 

Note  that  if  two  adjacent  grid  points  change  the  phase  at  the  same  time,  the 
above  process  is  still  applicable  except  that  additional  iteration  process  between  those 
adjacent  grid  points  is  needed  to  find  their  values  simultaneously. 

2.4.2     Interfacial  Conditions 

The  interface  shape  in  our  simulation  is  not  known  as  a  priori,  instead,  it  is 
determined  as  a  part  of  the  entire  solution  process.   At  each  time  step,  in  order  to 
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Figure  2.27:  Illustration  of  change  of  the  phase  on  some  grid  points  following  the 
interface  movement,  (a)  Location  of  the  interface  at  the  time  level  n.  The  shaded 
region  designates  the  vapor  phase.  At  this  instant,  the  grid  point  (i,  j)  lies  in  the 
vapor  phase,  (b)  Location  of  the  interface  at  the  time  level  n  +  1  after  moving.  The 
grid  point  (i,  j)  now  lies  in  the  liquid  phase.  Those  grid  points  which  undergo  the 
change  of  the  phase  are  marked  by  open  circles.  Also  shown  in  (b)  is  the  schematic 
of  finding  the  value  at  the  (i,  j)  by  drawing  a  normal  probe  passing  through  (»,  j), 
the  value  at  the  end  of  the  normal  probe  is  obtained  by  using  a  bilinear  interpolating 
function  over  points  numbered  1  to  4. 


determine  the  correct  interface  location  which  satisfies  all  interfacial  conditions  (1.23), 
(1.24),  (1.25)  and  (1.26),  an  iterative  process  for  updating  the  interface  location  is 
needed  as  described  in  section  2.4.3  and  2.4.4. 

The  interfacial  temperature  is  specified  using  Eq.  (1.26),  the  vapor  pressure  at 
the  interface  is  p^,,  then  the  liquid  pressure  at  the  interface  can  be  determined  by 
Eq.  (1.24). 

The  energy  conservation  condition  Eq.  (1.25)  is  used  to  estimate  the  interface 
velocity  relative  to  vapor  velocity  at  the  interface  each  time,  thus  it  is  naturally 
satisfied.  The  mass  conservation  condition  Eq.  (1.23)  serves  for  determining  the 
liquid  phase  velocity  at  the  interface  when  the  interface  velocity  and  vapor  phase 
velocity  at  the  interface  are  known.    Both  liquid  and  vapor  phase  velocities  at  the 
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interface  are  assigned  as  the  boundary  conditions  respectively  for  solving  the  flow 
fields  inside  and  outside  the  vapor  bubble. 

2.4.3     Iterative  Process  for  Determining  the  Interface  Shape 

In  this  section,  we  will  describe  the  process  for  determining  the  interface  shape 
using  normal  stress  balance  coupled  with  flow  solver  via  an  iterative  process.  This 
is  the  most  important  technique  and  difficult  numerical  procedure  in  this  work.  The 
method  depicted  in  the  following  is  used  to  determine  the  bubble  deformation  caused 
by  the  momentum  balance  with  no  phase  change  occurs.  The  technique  implemented 
in  the  present  research  is  the  one  similar  to  that  proposed  by  Ryskin  &  Leal  (1984) 
in  their  method  using  body-fitted  moving  grid.  Here  we  implemented  it  in  the  fixed 
grid  which  results  in  a  more  efficient  overall  solution  procedure. 

If  phase  change  is  to  be  considered  and  contributes,  in  conjunction  with  the  mo- 
mentum balance,  to  the  shape  change  of  the  bubble,  additional  algorithm  is  required 
and  developed  in  the  present  study.  They  are  given  in  the  next  section. 

Given  an  initial  interface  shape  and  flow  field,  material  properties  are  assigned 
in  the  liquid  and  vapor  phases,  the  flow  field  is  then  solved  until  a  converged  flow  field 
is  obtained  when  the  interface  is  fixed.  This  flow  field  serves  as  the  starting  point 
for  subsequent  moving  interface  problems.  Then  the  solution  algorithm  proceeds 
iteratively  through  the  following  steps  starting  from  iteration  number  k=0  for  next 
time  step  n  +  1,  the  initial  guess  for  flow  variables  at  time  step  n  +  1  is  <f)n+l>°  =  0n, 
the  initial  interface  velocity  is  w"^1'0  =  w"n(. 

The  interface  shape  is  determined  by  the  normal-stress  balance  Eq.  (1.24).  The 
technique  proposed  by  Ryskin  &  Leal  (1984)  is  adopted  to  obtain  the  estimate  of 
the  interface  shape  for  next  iteration.  The  idea  behind  this  technique  is  very  simple. 
This  approach  essentially  takes  the  local  imbalance  between  the  total  normal  stress 
rn,  which  includes  both  static  and  dynamic  pressure  and  viscous  contributions,  and 
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capillary  forces 

n(s)  =  rn-^-eK  (2.44) 

as  a  "driving  force"  which  causes  a  local  displacement  of  the  interface  in  the  normal 
direction,  the  magnitude  of  the  local  displacement  is  proportional  to  II (s).  Thus 

n+l,k+l  n+l,k 


Vint 


fc+1    =    yl+tl'k  +  pnk(s).ny  (2.45) 


where  the  under-relaxation  coefficient  j3  has  to  be  determined  by  numerical  experi- 
ment, its  typical  values  being  O(10~4  to  10-3)  in  our  computations. 

In  the  cases  where  interface  completely  encloses  one  of  the  two  fluids,  e.g.,  bubble 
in  our  study,  the  local  increment  in  the  location  of  interface  marker  points  must  be 
done  in  such  a  way  as  to  satisfy  the  mass  conservation  constraint,  i.e.  preserve  the 
bubble  volume  if  density  is  constant.  We  know  the  change  in  the  volume  between 
the  kth  and  (k  +  \)th  iterations  is 

f  IOC1**1  -  &*?  +  (ylV'k+1  -  y-:thknl2ds  (2.46) 

Jo 
where  the  integral  is  taken  along  the  interface.  Since 

K*£w  -  jfiW  +  (y%1M1  -  CW  -  nk(s)  (2.47) 

we  have  the  mass  constraint 

nk(s)ds  =  0  (2.48) 


f 

Jo 


10 

This  constraint  determines  a  free  constant  of  the  pressure  field  in  the  total  normal 
stress  rn  at  each  iteration  k  so  Eq.  (2.48)  is  satisfied. 

Even  after  this  constraint  has  been  satisfied,  the  bubble  may  still  change  volume 
slightly  at  each  iteration  due  to  higher-order  numerical  errors.  To  prevent  these 
small  changes  from  accumulating  and  becoming  significant,  a  simple  scaling  of  the 
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interface  can  be  done  following  the  Eq.  (2.45).  The  uniform  scaling  magnitude  A  of 
the  interface  location  in  the  normal  direction  can  be  determined  from 

/  Alds  =  AV  (2.49) 

Jo 

where  AV  is  the  error  in  the  volume.  This  simple  scaling  is  very  effective,  normally 

one  or  two  iterations  of  this  process  is  sufficient  to  bring  the  percentage  error  in  the 

volume  down  to  10"5. 

The  global  iterative  process  within  each  time  step  is  therefore  as  follows. 

1.  For  a  given  shape  of  the  bubble,  the  flow  field  is  computed  by  solving  the 
Navier-Stokes  equations  with  a  small  number  of  iterations  on  pressure  Poisson 
equation. 

2.  Knowing  the  flow  field,  the  normal  stress  balance  at  the  interface  is  calculated 
and  checked.  If  it  is  not  satisfied,  the  interface  shape  is  modified  according 
to  Eq.  (2.45)  so  as  to  reduce  the  imbalance  between  the  total  stress  and  the 
surface  tension  force. 

3.  After  each  interface  update,  the  mass  conservation  is  enforced  by  rescaling  the 
interface  using  Eq.  (2.49). 

4.  The  interface  normal  velocity  is  calculated  from  the  kinematic  condition  u™^1  = 
(x"^1  —  x"nt)  /At,  which  is  the  boundary  conditions  for  solving  the  flow  field. 

5.  Return  to  step  1  and  repeat  until  all  equations  and  boundary  conditions  are 
satisfied  to  a  predetermined  level  of  accuracy. 

2.4.4     Determining  Interface  Shape  with  Phase  Change 

The  interface  movement  when  phase  change  is  present  can  be  decomposed  into 
two  components:  local  propagation  velocity  due  to  evaporation  or  condensation,  and 
the  global  movement  resulted  from  buoyancy  when  there  is  no  phase  change. 
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The  interface  velocity  component  due  to  phase  change  effect  in  the  normal  di- 
rection, which  is  also  the  velocity  component  relative  to  the  vapor  phase,  is  as  Eq. 
(1.25),  i.e., 


\Un)int,p  ~  pe 


m    (kv\  (drv 


(2.50) 


dn       \ki  J  \  dn  t 

where  the  superscript  'n+1'  on  the  left  hand  side  denotes  the  n  +  1  time  step  while 

the  right  hand  side  uses  the  field  values  at  time  step  n.  The  subscript  'p'  on  the  left 

hand  side  means  this  is  the  interface  velocity  component  due  to  phase  change  effect. 

So  the  interface  movement  in  the  normal  direction  due  to  phase  change  is 

<&  •  »  =  <f  n  +  «£  •  n)  At  (2.51) 

After  we  have  this  new  location  of  the  interface,  we  can  use  the  process  in  the 
section  2.4.3  to  further  determine  the  shape  satisfying  the  momentum  balance  condi- 
tion. The  final  interface  velocity  at  time  step  n  +  1  obtained  from  kinematic  condition 
is  the  combination  of  the  components  owing  to  phase  change  and  momentum  balance. 

Hence,  the  vapor  phase  velocity  at  the  interface  in  the  normal  direction  is  given 
by: 

which  is  the  boundary  condition  for  solving  the  flow  field  in  the  vapor  bubble. 

The  liquid  phase  velocity  at  the  interface  in  the  normal  direction,  according  to 
Eq.  (1.23),  is 

Mr  =  Kffi1  *  [l  "  (  *)]  +  (£)  {un)Tl  (2.53) 

which  is  the  boundary  condition  for  liquid  phase. 

2.4.5     Discontinuity  of  Material  Properties 

In  the  currently  available  numerical  methods  reported  in  the  literature  for  sim- 
ulating two-phase  flow,  most  of  the  single  field  approaches  encounter  the  numerical 
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instability  caused  by  large  property  jumps  of  the  two  phases,  the  maximum  density 
ratio  for  phase  change  problems  which  is  reported  so  far  in  the  prior  publications,  to 
our  best  knowledge,  is  100.  However,  for  liquid-vapor  phase  phenomena,  the  density 
ratio  often  tops  1000.  In  a  typical  situation  of  water  at  atmospheric  pressure,  the 
density  ratio  of  liquid  and  vapor  is  roughly  1600.  For  some  air  bubble  problems  with 
no  phase  change,  we  saw  publications  dealing  with  density  ratio  as  large  as  1000  but 
with  a  finite  thickness  interface  treatment  to  prevent  instability. 

In  the  single  field  approaches,  for  any  cell  around  the  interface,  different  proper- 
ties will  appear  in  the  same  discretized  equation  on  that  cell  which  cause  the  stiffness 
in  that  equation  so  as  to  deteriorate  the  numerical  stability  although  easy  implemen- 
tation of  the  single  field  algorithm  is  quite  attractive. 

In  our  approach,  because  we  employ  the  explicit  interface  tracking  and  cut  cell 
technique  around  the  interface,  and  our  carefully  designed  computational  procedure 
solves  separately  the  two  fields  corresponding  to  two  phases  while  only  on  the  inter- 
face, the  two  fields  communicate  with  each  other  to  adjust  there  own  flow  field  to 
satisfy  the  interfacial  conditions,  the  fundamental  difference  from  other  approach  is 
that  the  discretized  equation  for  any  cell  including  those  around  the  interface  will 
involve  identical  material  properties  corresponding  to  the  phase  that  particular  cell 
lies  in.  Therefore,  there  is  no  stiffness  in  the  discretized  equation  on  any  cell,  the 
large  density  jump  will  not  be  expected  to  exert  severe  effect  on  the  numerical  stabil- 
ity in  our  approach  although  the  interface  has  zero  thickness.  This  benefit  is  at  the 
cost  of  complex  algorithm  design  and  coding  involved  in  the  development  of  present 
techniques. 

2.5     Discretization  on  Axisymmetric  Geometry 

In  the  simulations,  the  computational  domain  is  axisymmetric  where  it  is  as- 
sumed that  d/d9  =  0  and  the  swirl  velocity  ug  is  zero. 
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In  a  differential  formulation,  the  2-D  dimensional  conservation  equations,  Eqs. 
(1.4)-(1.7)  written  in  an  axisymmetric  cylindrical  coordinate  system,  can  be  further 
expressed  as 


1  d  /      a      1  d  ,      \ 

-—(rur)  +  -—(ruz)    =    0 

r  or  r  oz 
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or  equivalently  as 
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(2.55) 


If  the  coordinates  z  and  r  of  the  cylindrical  coordinate  system  are  replaced  by 
x  and  y,  the  analogy  between  Eq.  (2.54)  and  the  equations  in  Cartesian  coordinates 
becomes  apparent.  If  r  is  also  set  to  1,  Eq.  (2.54)  becomes  identical  to  those  in 
Cartesian  coordinates  except  with  an  additional  term  —/j,ur/r2  with  uz  —  ux  and 
ur  =  uy  for  differential  forms. 

In  terms  of  the  finite  volume  discretization  procedure,  the  conservation  equations 
written  in  integral  form  can  be  obtained  from  Eq.  (2.55)  by  integrating  the  both  sides 
over  an  axisymmetric  control  volume.  The  integral  equations  are  the  same  as  those  for 
planar  2D  cases  in  Cartesian  coordinate  except  for  an  additional  term  —pur/r2  ■  AV 
in  the  uT  momentum  equation,  where  AV  is  the  volume  of  the  control  volume.  Any 
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method  described  in  previous  chapters  are  equally  applicable  for  discretizing  the  2D 
axisymmetric  equations  in  cylindrical  coordinates.  However,  some  differences  need 
to  be  addressed  in  discretizing  the  equations  in  an  axisymmetric  control  volume. 

In  addition  to  the  surface  integrals  of  pressure  force  over  cell  faces  Z{  =  const, 
Zi-\  =  const,  Tj  —  const,  fy_i  =  const,  as  was  the  case  in  planar  2D  problems,  the 
radial  component  of  pressure  forces  onto  the  two  faces  9  =  const  has  to  be  considered 
also. 

The  other  difference  compared  to  planar  2D  problems  is  the  calculation  of  cell 

face  areas  and  cell  volume.   Since  the  control  volume  size  in  the  ^-direction  is  one 

radian,  the  areas  of  cell  faces  are  computed  in  the  same  way  as  in  the  planar  2D  case 

with  the  length  in  ^-direction  taken  into  account,  i.e.,  inclusion  of  a  factor  r/c  (where 

Tfc  denotes  the  cell  face  center) .  The  areas  of  the  front  and  back  faces  are  calculated 

in  the  same  way  as  the  volume  in  the  planar  geometry.  The  volume  of  axisymmetric 

CVs  with  any  number  of  faces  is  (see  Ferziger  &  Peric,  1996)  : 

,   nv 
AV  =  e  E  (*-i  -  t)  (r?-i  +  r*  +  r<r*-0  (2-56) 

where  Nv  denotes  the  number  of  vertices,  counted  counterclockwise,  with  i  =  0 
corresponding  to  j  =  Nv. 

With  the  similarity  of  the  governing  equations  for  the  planar  2D  problems  in 
Cartesian  coordinate  and  axisymmetric  2D  problems  in  cylindrical  coordinate,  the 
same  computer  code  can  be  used  for  both  planar  and  axisymmetric  2D  flows;  for  ax- 
isymmetric problems,  one  sets  r  —  y  to  account  for  the  length  in  the  third  dimension 
of  the  axisymmetric  control  volume,  for  planar  problems,  one  sets  r  =  1  because  the 
length  of  the  third  dimension  is  unity  for  planar  2D  control  volume. 


CHAPTER  3 
VALIDATION  OF  CARTESIAN  GRID  METHOD 

In  Chapter  2,  we  detailed  the  numerical  method  for  dealing  with  fixed,  moving, 
deforming  geometry  in  the  underlying  Cartesian  grid  system.  Also  the  phase  change 
effect  is  incorporated  into  the  moving  interface  updating  procedure.  The  moving 
interface  algorithm  is  based  upon  the  fundamental  Cartesian  grid  method  for  fixed 
complex  geometries  using  cut  cell  discretization.  Therefore,  validating  our  numerical 
method  is  a  two-step  process,  first,  the  spatial  discretization  scheme  using  Cartesian 
grid  method  for  fixed  geometries  is  to  be  checked.  Secondly,  on  top  of  the  fixed 
geometry  Cartesian  grid  method,  the  algorithm  for  moving  and  updating  dynamic 
interfaces  has  to  be  checked. 

It  is  the  task  of  this  chapter  to  fulfill  the  first  step  validation  of  the  overall 
numerical  method,  i.e.,  to  examine  the  accuracy  and  fidelity  of  the  fixed  interface 
Cartesian  grid  method  described  in  2.3.  In  the  next  Chapter,  we  will  make  some 
checks  on  the  overall  moving  interface  algorithm  involving  techniques  described  in 
2.4,  which  completes  our  second  step  validation  on  the  numerical  techniques. 

3.1     Validation  of  Global  Accuracy  of  Cut  Cell  Approach 

As  we  have  shown  in  section  2.3,  a  second-order  accurate  spatial  discretization 
scheme  adopting  cut  cell  approach  is  employed  in  the  vicinity  of  the  interface.  To 
confirm  that  we  obtain  the  nominal  second-order  accuracy  when  solving  the  flow  field 
over  complex  geometries,  the  most  straightforward  way  is  to  compute  a  flow  which 
has  a  curved  immersed  interface  and  one  for  which  an  analytical  solution  exists. 
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The  flow  chosen  here  corresponds  to  the  two-dimensional  Stokes  flow  past  a  circular 
cylinder  placed  next  to  a  moving  wall  (see  Figure  3.1). 

The  exact  solution  to  this  flow  was  given  by  Wannier  (1950)  and  is  reproduced 
here.  Here  we  have  simulated  this  flow  using  our  solver  on  four  different  uniform 
meshes.  The  meshes  have  equal  spacing  in  the  x  and  y  directions  and  have  Nx  and 
Ny  points  in  these  two  directions,  respectively.  The  global  error  in  the  numerical 
solution  is  computed  as 


N.N, 

—  y  i# 


NxNy  j=l 


-<t>7 


(3.1) 


In  order  to  simulate  Stokes  flow,  the  convection  terms  have  been  turned  off  in 
our  simulation.  Computations  have  been  carried  out  in  the  domain  shown  in  Figure 
3.1  with  the  exact  solution  imposed  on  the  boundaries.  In  Figure  3.2,  we  show  a 
log-log  plot  of  the  error  for  both  velocity  components  u  and  v  versus  Nx.  Also  shown 
is  a  line  with  a  slope  of  -2  which  corresponds  to  second-order  accurate  convergence. 


Figure  3.1:    Computational  domain  and  computed  streamline  pattern  of  Wannier 
flow. 
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Figure  3.2:  Global  error  of  u  and  v  as  a  function  of  mesh  points  for  Wannier  flow. 


The  plot  clearly  shows  that  the  global  error  in  our  computed  solution  decreases  in  a 
manner  consistent  with  a  second-order  accurate  scheme.  This  test  therefore  proves 
that  the  current  approach  of  treating  the  fluxes  in  the  boundary  cells  does  indeed 
result  in  a  solver  which  is  globally  second-order  accurate. 

3.2     Fidelity  of  The  Numerical  Methodology 

The  exact  solution  of  the  Wannier  flow  allows  us  to  confirm  the  accuracy  of  the 
solver  in  the  Stokes  flow  regime.  Here  we  validate  the  solver  in  the  finite  Reynolds 
number  regime  by  simulating  steady  and  unsteady  flow  past  a  circular  cylinder  im- 
mersed in  an  unbounded,  uniform  flow  over  a  wide  range  of  Reynolds  numbers  where 
the  Reynolds  number  is  defined  as  Red  =  V^djv  with  d  being  the  cylinder  diam- 
eter and  [Too  the  freestream  velocity.  This  flow  had  been  studied  quite  extensively 
in  the  past  and  a  number  of  numerical  and  experimental  datasets  exist  for  this  flow 
which  are  useful  for  the  purpose  of  validation.  Simulations  have  been  performed  at 
Red  =  20, 40, 80,  &  300  and  results  were  compared  with  established  experimental  and 
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Figure  3.3:  Typical  non-uniform  mesh  for  simulation  of  flow  past  a  circular  cylinder. 


numerical  results.  All  these  simulations  have  been  performed  in  a  large  30d  x  30d 
domain  so  as  to  minimize  the  effect  of  the  outer  boundary  on  the  development  of  the 
wake.  Figure  3.3  shows  the  152  x  156  non-uniform  mesh  used  in  the  low  Reynolds 
number  simulations.  At  the  inlet  and  along  the  top  and  bottom  boundaries  we  spec- 
ify velocities  corresponding  to  the  potential  flow  past  a  cylinder  and  a  homogeneous 
Neumann  boundary  condition  is  applied  at  the  exit  boundary.  We  have  also  tested 
larger  domain  sizes  in  order  to  ensure  that  the  results  presented  here  are  independent 
of  the  domain  size.  For  all  these  simulations  we  first  impose  a  small  asymmetric  dis- 
turbance at  the  inflow  boundary  for  a  short  period  of  time  and  then  allow  the  flow 
to  evolve  naturally  after  this.  For  Re^  =  20  and  40  the  wake  eventually  attains  a 
steady  symmetric  state  and  this  is  consistent  with  the  well  established  result  that 
the  cylinder  wake  is  stable  to  perturbations  below  Red  =  46  ±  1  (see  Jackson,  1987; 
Provensal  et  al.,  1987;  Williamson,  1996).  Once  the  flow  has  reached  a  steady  state 


95 


Figure  3.4:  Streamline  plot  of  flow  past  a  circular  cylinder,  (a)  Red  =  20  (b)  Red  =  40 


we  compute  the  drag  coefficient  defined  by  Co  =  nmiu^d  anc^  tne  length  of  the 
recirculation  zone  and  then  compare  these  with  established  results. 

The  streamline  plots  in  Figure  3.4a  and  b  show  the  mean  recirculation  regions 
behind  the  circular  cylinder  at  Red  ~  20  and  40  respectively.  In  Table  3.1,  results 
in  this  steady  flow  regime  obtained  by  the  current  method  are  compared  with  nu- 
merical simulation  by  Dennis  &  Chang  (1970)  as  well  as  experimental  measurements 
of  Tritton  (1959).  It  is  found  that  our  results  compare  well  with  other  numerical 
simulations  and  experiments. 

It  is  generally  accepted  that  the  wake  of  a  cylinder  immersed  in  a  freestream  first 
becomes  unstable  to  perturbations  at  a  critical  Reynolds  number  of  about  Red  = 
46  ±  1  (see  Jackson,  1987;  Provensal  et  ai,  1987).  Above  this  Reynolds  number, 
a  small  asymmetric  perturbation  in  the  near  wake  will  grow  in  time  and  lead  to 
an  unsteady  wake  and  Karman  vortex  shedding.    This  is  indeed  what  we  find  for 
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Table  3.1:  Comparison  of  mean  drag  coefficient,  length  of  wake  bubble  Lw  (measured 
from  rear  end  of  the  cylinder)  and  Strouhal  number  with  established  results. 


Reynolds  Number  —*■ 

20 

40 

80 

300 

Mesh  — > 

152  x  156 

152  x  156 

217  x  183 

217  x  183 

Study  4, 

cD 

Lw/d 

cD 

Lw/d 

cD 

St 

cD 

St 

Tritton  (1959) 

2.22 

— 

1.48 

— 

1.29 

Weiselsberger  (1922) 

2.05 

— 

1.70 

— 

1.45 

— 

1.22 

— 

Dennis  &  Chang  (1970) 

2.05 

0.94 

1.52 

2.35 

Fornberg  (1980) 

2.00 

0.91 

1.50 

2.24 

Williamson  (1996) 

0.15 

— 

0.20 

Current 

2.03 

0.92 

1.52 

2.27 

1.37 

0.15 

1.38 

0.21 

our  simulation  at  Red  —  80.  Figure  3.5  shows  the  variation  of  the  lift  and  drag 
coefficients  with  time  and  it  shows  how  vortex  shedding  develops  to  a  periodic  state 
in  time.  The  computed  mean  drag  coefficient  from  the  current  simulation  is  about 
1.37  which  lies  between  the  two  experiments  results.  The  Strouhal  number  for  vortex 
shedding  is  defined  as  St  =  fd/Uoo,  where  /  is  the  shedding  frequency  and  is  one 
of  the  key  quantities  that  characterizes  the  vortex  shedding  process.  Here  we  have 
estimated  the  Strouhal  number  from  the  periodic  variation  of  the  lift  coefficient  and 
the  value  comes  out  to  be  0.15  which  compares  very  well  with  the  value  obtained 
from  experiments  (see  Williamson,  1996).  Figure  3.6  shows  a  plot  of  the  streamline 
pattern  and  spanwise  vorticity  contour  at  one  instant  and  both  plots  show  the  classical 
Karman  vortex  street. 

In  addition  to  the  low  Reynolds  number  simulations,  we  have  also  carried  out 
a  simulation  at  a  moderately  high  Reynolds  number  of  300.  This  simulation  serves 
to  demonstrate  that  the  current  methodology  is  capable  of  resolving  thin  boundary 
layers  that  develop  in  flows  at  these  Reynolds  numbers.  The  mesh  used  for  this 
simulations  is  the  same  as  that  used  for  Red  =  80.  This  is  a  relatively  coarse  mesh  and 
this  coarse  resolution  severely  tests  the  discretization  scheme  used  in  our  solver  for  the 
boundary  cells.  The  mean  drag  and  Strouhal  numbers  have  been  computed  from  time 
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Figure  3.5:  Variation  of  lift  and  drag  coefficients  with  time  for  Red  =  80 


Figure  3.6:  Instantaneous  streamline  plot  and  vorticity  contour  plot  in  the  near  wake 
of  the  circular  cylinder  for  Red  =  80  (a)  Streamline  plot  (b)  Vorticity  contours 
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variation  of  lift  and  drag  coefficients  and  they  are  included  in  Table  3.1.  It  can  be  seen 
that  the  Strouhal  number  matches  well  with  the  experiments  of  Williamson  (1996).  It 
should  be  pointed  out  that  at  this  Reynolds  number  the  cylinder  wake  is  intrinsically 
three-dimensional  whereas  our  simulation  is  two-dimensional  and  therefore  does  not 
allow  spanwise  variations.  It  is  well  known  that  the  2-D  simulation  of  a  flow  which 
is  really  3-D  over-predicts  the  drag.  This  is  indeed  what  we  observe  for  the  current 
simulations.  The  drag  coefficient  from  our  2-D  simulations  is  about  12%  higher  than 
the  experimentally  determined  value  of  Wieselsberger  (1922). 

In  Figure  3.7,  we  have  shown  contour  plots  of  spanwise  vorticity  at  one  instant. 
Figure  3.7a  shows  a  view  of  of  the  wake  that  extends  to  about  10  diameters  down- 
stream from  the  cylinder  and  as  expected,  this  plot  shows  the  formation  and  evolution 
of  compact  Karman  vortices  in  the  wake.  Figure  3.7b  is  a  closeup  view  of  the  flow 
around  the  cylinder  and  the  mesh  superposed  on  the  greyscale  contour  plot  clearly 
shows  that  there  are  fewer  than  five  points  in  the  attached  boundary  layer.  It  can 
be  seen  that  even  with  a  relatively  low  resolution  provided  here,  the  boundary  lay- 
ers on  the  cylinder  surface  are  smooth  indicating  that  the  current  treatment  of  the 
boundary  cells  adequately  resolves  thin  boundary  layers. 

3.3      Summary 

In  this  chapter,  the  numerical  algorithm  employing  Cartesian  grid  method  to 
simulate  flow  with  fixed  complex  geometries  are  validated  through  extensive  test 
cases.  The  results  have  confirmed  the  accuracy  and  fidelity  of  the  current  compu- 
tational method.  The  first  step  in  our  validation  process  is  completed.  In  the  next 
chapter,  we  will  conduct  the  second  step  of  the  validation  process,  i.e.,  validate  nu- 
merical techniques  for  moving  interfaces. 
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Figure  3.7:  Spanwise  vorticity  contour  plots  in  the  wake  of  the  circular  cylinder  for 
Rej,  =  300  (a)  View  extending  to  9d  downstream  of  the  cylinder,  (b)  Closeup  view 
showing  the  resolution  provided  to  the  attached  boundary  layers  and  separated  shear 
layers 


CHAPTER  4 
BUOYANCY-DRIVEN  BUBBLE  MOTION 

The  three  major  components  of  the  present  numerical  methodology  are  interface 
representation  scheme,  Cartesian  grid  flow  solver  and  moving  interface  algorithms. 
In  previous  chapters,  the  capability  of  two  of  the  three  components  are  examined. 
The  accuracy  of  the  interface  tracking  algorithms  is  tested  and  demonstrated.  The 
accuracy  and  fidelity  of  the  Cartesian  grid  method  for  solving  fluid  flows  over  fixed, 
sharp  interfaces  with  complex  geometries  are  also  validated  through  extensive  cases. 

Since  the  main  objective  of  this  research  is  to  study  vapor  bubble  heat  transfer 
by  direct  numerical  simulations,  we  present  in  this  chapter  the  development  of  pre- 
vious numerical  methods  on  an  axisymmetric  cylindrical  coordinates  system.  All  the 
vapor  bubble  simulations  are  performed  on  this  system.  Specifically,  the  solution  pro- 
cedure dealing  with  the  axisymmetric  geometry,  the  B-spline  curve  fitting  and  fairing 
algorithms  for  interface  representation  and  curvature  evaluation,  the  interface  up- 
dating procedure  based  on  normal  stress  balance,  and  the  global  iterative  process  for 
solving  flow  fields  and  interface  motion  simultaneously  are  validated.  The  problems 
simulated  in  this  chapter  involve  no  phase  change.  The  flow  investigated  is  unsteady, 
axisymmetric  and  laminar.  For  this  type  of  problems,  there  have  been  many  reported 
investigations.  The  results  abound  and  physical  aspects  are  well  understood  so  that 
we  may  validate  the  entire  numerical  method  by  comparing  our  solutions  with  others. 

On  the  other  hand,  the  buoyancy-driven  motion  of  a  bubble  is  a  typical  moving 
boundary  problem  in  fluid  mechanics.  The  physical  process  is  highly  complex  subject 
to  multiple  effects.  By  simulating  such  flows  with  the  dynamic  interface,  we  can 
obtain  a  valuable  understanding  and  insight  of  the  solution  characteristics  -  both  the 
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flow  field  and  bubble  deformation  -  of  our  method.  These  simulations  thus  serve  as 
the  essential  precursor  to  numerical  simulations  of  rising  bubble  with  phase  change 
that  will  be  detailed  in  the  following  chapters. 

4.1     Summary  of  Previous  Study 

Although  a  single  bubble  rising  in  a  liquid  has  been  studied  for  a  number  of 
years,  most  theoretical  results  are  limited  to  very  small  deformation  at  either  low  or 
high  Reynolds  numbers.  For  example,  at  very  low  Reynolds  numbers,  there  exists 
the  theoretical  model  by  Taylor  &  Acrivos  (1964)  based  on  an  asymptotic  theory. 
At  high  Reynolds  numbers,  only  boundary-layer  approximations  (see  Moore,  1963; 
Harper  &  Moore,  1968)  and  semi-empirical  models  (see  Parlange,  1969)  are  available. 
All  of  the  above  investigations  assume  that  the  bubble  maintains  a  spherical  shape 
which  is  rather  unrealistic  at  high  Reynolds  or  Weber  numbers. 

Experimental  investigation  of  Saffman  (1956)  and  Bhaga  &;  Weber  (1981)  pro- 
vided a  fairly  detailed  picture  of  the  gas  bubble  upward  motions  in  a  liquid. 

Ryskin  &  Leal  (1984)  have  reported  the  first  successful  theoretical  solution  for 
motion  of  bubbles  with  finite  degree  of  deformation  using  body-fitted,  moving  grid 
techniques.  In  their  model,  the  interface  is  treated  as  zero  thickness;  the  shape  is 
explicitly  determined  by  the  stress  balance  at  the  interface.  However,  the  problems 
they  considered  involve  one  viscous  fluid  surrounding  the  bubble  and  a  void  bubble 
with  pv  =  0,  fiv  =  0.  There  is  no  flow  field  inside  the  bubble.  So  the  interfacial 
conditions  involve  forces  on  only  one  side  of  the  bubble  interface.  Later,  Dandy  & 
Leal  (1989)  have  extended  the  numerical  method  of  Ryskin  &  Leal  (1984)  to  consider 
the  deformable  drop  problems  involving  two  viscous  fluids  both  inside  and  outside 
of  the  drop.  In  both  Ryskin  &  Leal  (1984)  and  Dandy  &  Leal  (1989),  only  steady- 
state  problems  are  considered.  Furthermore,  orthogonal,  body-fitted  coordinates  are 
adopted  to  generate  the  grid  system. 
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In  recent  years,  many  numerical  simulations  of  the  bubble  motion  have  been 
reported  in  the  literature.  However,  in  virtually  all  cases,  the  interface  is  not  sharply 
defined,  and  the  stress  balance  is  enforced  via  several  cell  spacings,  instead  of  at 
precisely  defined  locations.  Examples  include  the  immersed  boundary  method  ,and 
volume  of  fluids  method  as  listed  in  Chapter  1.  The  nature  of  those  methods  does  not 
render  possibly  the  accurate  prediction  of  the  exact  location  of  deformed  interface. 
While  our  method  solves  the  two  separate  fields  with  stress  conditions  matched  at 
the  interface  to  obtain  the  exact,  instantaneous  interface  shape.  The  comparison 
between  our  results  and  those  obtained  recently  (e.g.,  see  Chen,  Garimella,  Reizes  & 
Leonardi,  1999)  is  therefore  not  meaningful  owing  to  the  difference  in  dealing  with 
interface  stress  conditions. 

To  help  assess  the  performance  of  the  present  method,  we  will  make  direct  com- 
parison with  the  results  reported  by  Ryskin  &  Leal  (1984)  for  the  bubble  motion  at 
different  Reynolds  and  Weber  numbers.  It  is  noted  that  the  results  of  Ryskin  &;  Leal 
(1984)  agree  well  with  experimental  study  of  Saffman  (1956)  and  Bhaga  &  Weber 
(1981). 

4.2     Grid  Resolution  Study  and  Convergence  Criteria 

Referring  to  Figure  4.1,  unless  otherwise  mentioned,  all  the  cases  reported  in  the 
above  are  computed  on  a  10  x  3  domain  with  a  250  x  75  uniform  mesh.  The  boundary 
conditions  for  the  cases  presented  below  are:  at  three  far  sides  of  the  boundary,  the 
outflow  (zero  velocity  gradient)  condition  is  specified  for  velocity.  At  the  line  of 
symmetry,  the  mirror  condition  for  all  variables  is  used. 

The  effect  of  grid  resolution  on  solution  accuracy  is  examined  first.  Unless 
otherwise  mentioned,  all  simulations  in  this  study  employ  the  same  grid  resolution 
around  the  interface,  that  is,  the  number  of  cells  across  the  initial  bubble  diameter 
is  25.  In  order  to  assess  the  grid  dependency  of  the  solution,  we  have  conducted 
computations  for  one  case:  Re  =  100,  We  =  4,  pv/pi  =  0.001,  fa/tk  =  10,  using 
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Figure  4.1:  Schematic  of  the  computational  domain  and  boundary  conditions.  For 
cases  involving  stationary  bubble,  L/R0  =  48,  H/Rq  —  24;  For  rising  bubble  cases 
without  phase  change,  L/Rq  =  20,  H/R0  —  6;  For  rising  bubble  cases  with  phase 
change,  L/Rq  =  20,  H/Rq  =  8 

500  x  150  grid  with  50  cells  across  the  initial  bubble  diameter.  Figure  4.2  shows  the 
time  history  of  the  aspect  ratio  of  the  bubble.  The  aspect  ratio  is  defined  as  the 
length  along  the  major  axis  divided  by  that  along  the  minor  axis.  As  can  be  seen  in 
Figure  4.2,  the  difference  between  the  results  from  the  two  grids  is  small. 

The  criteria  for  determining  the  convergence  are:  the  maximum  norm  of  the 
absolute  error  must  be  less  than  10-6  for  momentum  equations  (2.12);  10-4  for  pres- 
sure Poisson  equation  (2.15);  10-6  for  energy  equation  in  phase  change  cases;  10~3 
for  evaluating  the  normal  stress  balance  (2.45);  When  these  convergence  criteria  are 
satisfied,  the  maximum  norm  of  the  residual  of  mass  conservation  for  the  entire  flow 
field  is  less  than  10"3. 

4.3     Results  For  Buoyancy  Driven  Bubble  Motion 

Representative  cases  of  Ryskin  &  Leal  (1984)  are  simulated  in  the  present  study. 
We  compute  the  solutions  for  Re  in  the  range  of  1  <  Re  <  100,  and  We  from  0 
up  to  20  for  Re  <  20  and  up  to  10  for  Re  >  50.  The  parameter  ranges  are  in  line 
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Figure  4.2:   Comparison  of  the  aspect  ratio  history  for  the  rising  bubble  for  Re  = 
100,  We  =  4,  pv/pi  =  0.001,  pv/pi  =  1.0  on  the  250  x  75  grid  and  the  500  x  150  grid 


with  those  of  Ryskin  &  Leal  (1984).  All  computations  are  done  in  time  dependent 
manner. 

To  facilitate  the  direct  comparison,  the  condition  of  balance  between  the  drag 
force  and  buoyancy  force  used  in  the  work  of  Ryskin  &  Leal  (1984)  is  employed  in 
the  present  study  to  determine  the  Froude  number  for  the  hydrostatic  pressure  term 
of  (1.24).  The  relation  used  by  Ryskin  &  Leal  (1984)  is 


2Rg  _  3 

TP~4Cd 


(4.1) 


where  R  is  the  bubble  radius,  U  is  the  terminal  velocity  of  the  bubble,  Co  is  the 
drag  coefficient.  The  left  hand  side  of  (4.1)  is  actually  1/Fr.  For  each  case,  we  use 
(4.1)  to  find  the  Froude  number  from  a  given  Cd  value.  This  procedure  ensures  that 
the  scaling  processes  between  the  current  and  that  used  in  Ryskin  &  Leal  (1984)  are 
consistent.  Of  course,  the  drag  coefficients  are  computed  from  the  solution  obtained. 
The  drag  coefficients  obtained  in  the  present  study,  given  in  Table  4.1,  are  those  when 
the  bubble  reaches  constant  rising  velocity. 
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The  drag  coefficients  obtained  by  our  simulations  and  by  Ryskin  &  Leal  (1984) 
are  summarized  in  Table  4.1.  Also  included  in  the  table  is  the  error  estimate  reported 
in  Ryskin  &  Leal  (1984),  based  on  energy  dissipation  analysis  of  their  numerical 
simulation.  This  information  helps  one  gain  a  sense  of  the  accuracy  level  in  that 
work. 

Table  4.1:  Comparison  of  drag  coefficients  from  present  simulations  (first  row) 
with  those  obtained  by  Ryskin  &  Leal  (second  row)  using  integration  of  the  forces 
at  the  surface.  The  third  row  shows  the  relative  deviation  of  drag  coefficients 
computed  via  energy  dissipation  in  Ryskin  &  Leal's  computations.  In  Ryskin  & 
Leal's  computations,  the  bubble  is  consider  to  be  a  void,  while  in  the  present  case, 
pt/pv  =  1605,  in/fly  =  22. 
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The  favorable  overall  agreement  in  drag  coefficients  between  the  two  simulations 
shows  that  our  method  is  capable  of  correctly  predicting  the  dynamic  behavior  of  a 
coupled  system  involving  the  liquid  flow  field  and  vapor  bubble. 
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The  computed  bubble  shapes  for  selected  cases  in  Table  4.1  are  shown  in  Figure 
4.3.  The  shapes  are  the  ones  when  the  unsteady  bubble  motion  reaches  the  terminal 
velocity. 
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Overall,  the  trend  of  bubble  shapes  changing  with  increased  We  is  in  agreement 
with  common  experimental  observation:  spherical  to  oblate-ellipsoidal  and  then  to 
oblate-ellipsoidal/spherical  cap  Bhaga  &  Weber  (1981). 

Figure  4.4  shows  development  of  the  bubble  shapes  for:  (a)  Re  —  5,  We  — 
10,  pi/pv  =  1605,  pi/pv  =  22,  (b)  Re  =  2,  We  =  16,  Pl/pv  =  1605,  pt/pv  =  22. 
The  corresponding  streamlines  for  the  two  cases,  plotted  based  on  the  coordinate 
fixed  at  the  middle  of  the  lower  surface  of  the  moving  bubble,  are  shown  in  Figure 
4.5.  Two  recirculating  structures  are  observed  in  each  case,  one  inside  the  bubble, 
and  the  other  caused  by  the  interaction  between  the  bubble  and  the  liquid.  It  is 
interesting  to  observe  with  Re  =  2  and  We  —  16,  the  recirculating  flow  in  the  liquid 
phase  is,  as  expected,  attached  to  the  bubble;  for  Re  =  5  and  We  ~  10,  it  tends  to 
detach  from  the  bubble. 

Figure  4.6  compares  the  steady  bubble  shapes  at  three  density  ratios.  The 
differences  observed  are  small.  The  reason  for  this  phenomenon  is  that  by  fixing  Re 
and  We,  the  only  impact  from  the  density  ratio  is  via  the  unsteady  and  convection 
terms  in  the  momentum  equation  in  the  vapor  domain.  For  a  rising  bubble,  since  the 
fluid  dynamics  inside  the  bubble  is  induced  by  the  interface  movement,  for  the  present 
Re  and  We  (defined  based  on  the  properties  of  the  liquid  phase),  the  impact  from  the 
vapor  portion  of  fluid  dynamics  is  limited.  Accordingly,  only  minor  differences  are 
observed  in  Figure  4.6.  Nevertheless,  it  is  reasonable  to  observe  that  as  the  density 
ratio  increases,  the  bubble  becomes  less  deformed.  Figure  4.7  shows  the  bubble 
shapes  at  selected  time  instants  for  the  three  density  ratios.  These  observations  have 
been  reported  previously  in  other  studies,  including  Dandy  &  Leal  (1989).  The  drag 
coefficients  for  these  three  cases  from  our  simulations  are  1.29,  1.32,  1.34  for  density 
ratios  of  0.1,  0.01,  0.001,  respectively.  The  drag  coefficient  reported  in  Dandy  k 
Leal  (1989)  for  varying  density  ratios  under  the  corresponding  Reynolds  and  Weber 
numbers  is  1.29. 
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Figure  4.4:    Development  of  bubble  shapes:    (a)  Re   =   2,    We   =   16,    pt/pv 
1605,  pi/pv  -  22,  and  (b)  Re  =  5,  We  =  10,  pt/pv  =  1605,  /x,/^,  =  22  . 
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Figure  4.5:  Flow  structures  at  the  terminal  state  for  cases  corresponding  to  Figure 
4.4,  (a)  Re  =  2,  We  =  16,  pt/pv  =  1605,  pt/pv  =  22,  and  (b)  Re  =  5,  We  =  10, 
Pi/pv  =  1605,  pi/pv  =  22.  The  streamlines  are  observed  on  the  reference  frame 
attached  to  the  moving  bubble. 
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Figure  4.6:   The  terminal  state  shapes  for  cases  with  Re  =  100,   We  —  4,   Fr 
1,  pv/pi  =  0.1,  0.01,  0.001,  pv/pt  =  1 


To  further  illustrate  the  effect  of  density  ratio  on  the  computational  performance, 
Figure  4.8  compares  the  convergence  histories  between  two  cases  with  different  density 
ratios.  The  residues  of  both  the  Young-Laplace  equation,  (1.24),  and  the  pressure 
Poisson  equation,  at  a  given  time  step,  are  shown.  The  residues  are  based  on  the 
sum  of  the  absolute  value  of  the  residue  computed  in  each  cell.  The  levels  shown  in 
Figure  4.8  are  not  normalized.  The  figure  demonstrates  that  the  present  method  is 
robust  in  terms  of  handling  the  large  property  variations  across  the  phase  boundary. 

Figure  4.9  shows  the  flow  structures  corresponding  to  three  density  ratios,  each 
with  Re  =  100,  We  =  4,  pv/pi  =  1.0.  This  figure  corresponds  to  the  same  pa- 
rameters as  those  shown  in  Figure  4.6  and  4.7.  In  all  cases,  the  recirculating  wake 
is  detached  from  the  bubble.  Again,  there  is  no  significant  difference  for  different 
density  ratios. 
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Figure  4.7:  The  shape  evolution  for  cases  with  Re  =  100,  We  =  4,  Fr  =  1,  fi*/m  — 
1  at  equal  time  intervals. 
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Figure  4.8:  Convergence  paths  of  the  Young-Laplace  equation  at  the  interface  and 
the  pressure  equation  in  the  entire  domain  within  a  time  step.  Here  two  different 
density  ratio  cases  are  shown  with  Re  =  100,  We  =  4,  Fr  —  1,  //„////  =  1. 
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Figure  4.9:  Flow  structures  at  the  terminal  state  for  cases  with  Re  =  100,  We  = 
4,  Fr  =  1,  fj,v/ni  =  1,  pv/pi  —  0.1,  0.01,  0.001;  The  streamlines  are  observed  on  the 
reference  frame  attached  to  the  moving  bubble. 
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4.4     Summary 

By  applying  the  present  method  to  simulating  the  buoyancy-driven  motion  of  a 
bubble,  the  method  is  proved  to  be  accurate  and  capable  of  resolving  the  interfacial 
dynamics  and  shapes.  It  should  be  noted  that  for  the  bubble  motion  considered 
here,  to  my  best  knowledge,  comparison  of  drag  coefficients  with  experiments  and 
established  theoretical  results  were  not  reported  in  other  fixed  grid  based  methods 
because  those  methods  usually  are  unable  to  resolve  the  exact  sharp  interface.  While 
the  drag  coefficient  is  strongly  influenced  by  interfacial  dynamics,  correct  calculation 
of  drag  coefficients  requires  the  method  to  resolve  the  dynamic  interface.  To  this 
extent,  the  present  method  is  suitable  for  problems  when  interfacial  dynamics  plays 
an  important  role  in  the  system. 


CHAPTER  5 

BUBBLE  RISE  AND 
GROWTH  DUE  TO  EVAPORATION 

In  Chapter  4,  the  simulation  results  based  on  the  current  method  for  buoyancy- 
driven  bubble  motion  without  phase  change  are  compared  with  other's  work.  In  this 
chapter,  the  results  for  the  heat  transfer-controlled  bubble  growth  due  to  evaporation 
in  a  superheated  liquid  under  either  zero  gravity  or  normal  gravity  conditions,  are 
presented. 

5.1     Stationary  Bubble  Growth 

Under  the  zero  gravity  condition,  the  bubble  motion  resembles  the  ideal  case 
of  a  stationary  bubble  growth  studied  in  many  early  works  such  as  Forster  &  Zu- 
ber  (1955),  Scriven  (1959)  and  Prosperetti  &  Plesset  (1978).  In  those  studies,  the 
bubble  is  assumed  in  a  spherical  shape,  thus  a  1-D  problem  with  the  bubble  radius 
as  the  dependent  variable  is  solved  in  conjunction  with  the  thermal  boundary  layer 
approximation.  The  momentum  effect  in  the  liquid  and  vapor  phases  on  the  bubble 
growth  and  its  shape  is  totally  discarded.  In  doing  so,  the  bubble  growth  rate,  i.e., 
the  radial  velocity  depends  only  on  the  Jakob  number.  All  those  studies  concluded 
that  the  time  evolution  of  the  bubble  growth  radius  follows  the  t1/2  law. 

For  the  heat  transfer-controlled  stationary  bubble  growth,  the  major  transport 
mechanism  is  the  heat  diffusion,  the  appropriate  scaling  for  the  velocity  will  be  based 
on  the  diffusion  mechanism:  ur  =  ai/L,  where  ol\  is  the  liquid  thermal  diffusivity, 
and  L  is  the  initial  bubble  diameter.  Using  the  diffusion  scale  for  the  velocity,  the 
Peclet  number  is  always  1.0. 
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5.1.1     Case  with  Thermal  Properties  of  Water 

First  we  computed  a  stationary  bubble  growth  based  on  the  thermal  properties 
of  water  under  the  atmospheric  conditions,  as  listed  in  Table  1.3.  A  AT  =  \°C 
superheat  is  considered.  Under  this  superheat,  the  bubble  departure  diameter  is 
about  0.3  mm  according  to  empirical  calculations  (see  Carey,  1992),  which  is  our 
initial  bubble  diameter,  and  thereby  the  length  scale  L. 

Using  the  thermal  properties  of  the  water  listed  in  Table  1.3,  the  thermal  diffu- 
sivity  of  water  at  atmospheric  pressure  is  a.\  =  1.679  x  10-7  m2/s.  Hence  the  velocity 
scale  is 

uT  =  at/L  =  (1.679  x  10"7)/(0.3  x  10"3)  =  5.  5967  x  10"4     m/s 

From  the  velocity  and  length  scales,  all  required  dimensionless  parameters  for 
this  case  are  estimated  as 

z?p  _  (HUrL  __   (958.3)(5.5967xl0-4)(0.3xl0-3) 

Ke  ~      w      -  277.53x10-6 -  0.  58 

Wr  _  flggl   _   (958.3)(5.5967xl0-4)2(0.3xl0-3)   _  6 

vv  e  —       a       —  0.0589  l,ugAW 

Jn  -  £LS££L  -   058.3)  (4220)(1)   _  3   Q 
0  u         pv       X  (0.597)2256700  °'  u 

The  Prandtl  number  is  1.72.  The  Weber  number  at  this  magnitude  essentially 
means  no  driving  force  for  the  deformation  because  everything  involved  is  radially 
symmetric. 

The  above  dimensionless  parameters  are  thereby  used  in  the  simulation.  At 
time  t  =  0,  the  liquid  phase  is  uniformly  superheated,  i.e.,  an  initial  temperature 
step  exists:  T(\x\  >  R0)  =  T^,  T(\x\  <  Rq)  =  Tv.  A  non-uniform  mesh  on  a  domain 
of  48i?0  x  24i?0  is  used  in  the  computations.  Rq  is  the  initial  bubble  radius. 

The  resulted  transient  dimensionless  bubble  radius  R(t)/Ro,  as  a  function  of  time 
is  plotted  in  Figure  5.1(a).  The  stationary  bubble  growth  exhibits  an  asymptotic  slope 
of  1/2  in  the  log-log  plot.  This  agrees  with  the  well  accepted  theory  that  the  growth 
radius  is  proportional  to  the  square  root  of  the  time.  The  1/2  law  is  valid  when  the 
thermal  boundary  layer  around  the  bubble  is  well  developed  after  the  initial  stage. 
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Figure  5.1:  Growth  rate  of  the  bubble  radius  for  a  stationary  bubble  obtained  in  the 
present  simulation  for  Re  =  0.58,  We  =  1.53  x  10-6,  Ja  =  3.0,  Pr  =  1.72,  Pe  =  1.0, 
Pi/Pv  =  1605,  [iiliiv  =  22.  (a)  Dimensionless  bubble  radius  (b)  Growth  portion  of 
the  bubble  radius. 
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In  Figure  5.1(b),  we  plotted  only  the  growth  portion  (R/R0  -  1.0)  of  the  bubble 
radius  versus  time.  The  plot  also  shows  that  the  bubble  growth  rate  approaches  the 
diffusion  controlled  limit  t1/2. 

The  above  case  with  water  required  approximately  20, 000  time  steps  for  the  bub- 
ble growth  rate  to  approach  the  diffusion  controlled  limit  because  of  the  magnitude 
of  the  Jakob  number  J  a  —  3.0,  which  requires  very  small  time  steps. 

5.1.2     Case  with  Thermal  Properties  of  Ethanol 

In  the  second  case,  we  use  the  thermal  properties  of  Ethanol  as  listed  in  Table 
1.4.  To  be  efficient,  a  smaller  Jakob  number  Ja  =  0.1  is  chosen  so  that  a  larger  time 
step  could  be  used.  All  the  other  dimensionless  parameters  are  determined  in  a  way 
similar  to  that  of  the  water  case  described  earlier.  The  Prandtl  number  is  Pr  =  8.37. 
All  dimensionless  parameters  are  listed  as 

Re  =  0.12 

Pe  =  1.0 

We  =  3.91  x  10-7 

Pr  =  8.37 

Ja  =  0.1 

Figure  5.2  shows  the  growth  portion  of  the  bubble  radius  versus  time  for  this 
case.  The  growth  radius  is  also  following  the  relation  R(t)  oc  t1/2. 

The  two  cases  confirm  the  analytical  prediction  for  the  growth  of  the  stationary 
bubble,  namely,  R(t)  oc  t1/2. 

5.2     Translating  Bubble  Growth 

Under  the  influence  of  gravity,  the  bubble  would  be  rising  and  growing  simulta- 
neously. Darby  (1964)  and  Ruckenstein  &  Davis  (1971)  suggest  that  the  growth  rate 
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Figure  5.2:  The  dimensionless  growth  radius  for  a  stationary  bubble  obtained  in  the 
present  simulation  versus  time  for  Re  —  0.12,  We  =  3.91  x  10-7,  J  a  =  0.1,  Pr  =  8.37, 
Pe  =  1.0,  pi/pv  =  527,  fj.i/fiv  =  41. 


is  significantly  higher  when  the  relative  motion  between  the  bubble  and  the  surround- 
ing liquid  is  present  based  on  their  analysis.  A  growth  rate  R(t)  oc  t2^  is  speculated 
under  such  conditions. 

For  the  rising  and  growing  bubble  due  to  buoyancy,  an  appropriate  scale  for  the 
velocity  is  uT  =  \fgL,  where  g  is  gravity  and  L  is  the  initial  bubble  diameter.  Using 
this  scale,  the  Froude  number  is  always  1.0. 

Two  cases  were  computed  for  the  rising  bubble  case  based  on  the  thermal  prop- 
erties of  water  and  ethanol,  respectively.  The  dimensionless  parameters  for  the  two 
cases  are  summarized  in  the  following.  The  simulations  were  conducted  in  the  time 
dependent  manner  as  in  all  other  cases  in  this  work.  The  simulations  were  stopped 
when  the  thermal  boundary  layer  around  the  bubble  became  fully  developed.  The 
dimensionless  times  when  the  two  cases  were  stopped  respectively  were  not  the  same. 
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1.  The  first  case  we  simulated  is  based  on  the  thermal  properties  of  water  under 
atmospheric  conditions.  However  for  convenience,  the  Prandtl  number  is  kept 
as  1.0.  Reynolds  number  is  10,  Peclet  number  is  10  because  of  the  value  of 
Prandtl  number  Pr  =  1.0.  The  Weber  number  is  set  as  0.2  which  corresponds  to 
a  slight  deformation  case.  The  Jakob  number  is  1.0.  The  system  dimensionless 
parameters  are  listed  in  the  following. 

Re  =  10.0 

Fr  =  1.0 

Pe  =  10.0 

We  =  0.2 

Pr  =  1.0 

Ja  =  1.0 

2.  The  second  case  we  simulated  is  based  on  the  thermal  properties  of  Ethanol 
under  the  atmospheric  conditions,  as  listed  in  Table  1.4.  A  AT  =  10°C  su- 
perheat is  considered.  Under  this  superheat,  the  bubble  departure  diameter 
is  about  1.0  mm  according  to  empirical  calculations  (see  Carey,  1992),  which 
is  our  initial  bubble  diameter,  and  thereby  the  length  scale  L.  The  velocity 
scale  is  thus  ur  =  yfgL  =  9.9  x  10~2  m/s.  With  these  reference  scales,  the 
dimensionless  parameters  for  this  case  are  estimated  using  the  procedure  given 
in  Section  5.1.1  and  listed  as  follows: 

Re  =  174.81 

Fr  =  1.0 

Pe  =  1463.7 

We  =  0.42 

Pr  =  8.37 
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Ja  =  16.43 

The  bubble  growth  radius  history  for  these  two  cases  is  shown  in  Figure  5.3. 
The  dotted  lines  have  slopes  1/2  and  2/3  respectively.  The  slopes  of  the  bubble 
growth  radius  from  the  present  simulation  fall  in  between  the  two  dotted  lines.  This 
figure  shows  that  the  growth  rate  is  higher  when  bubble  is  rising  as  compared  to  a 
stationary  bubble.  However,  the  growth  rate  is  approaching  but  does  not  reach  the 
t2/3  limit,  which  is  probably  due  to  the  deformation  of  the  bubble  from  a  sphere. 

At  the  end  of  the  simulations,  the  streamlines  for  the  two  cases,  plotted  based 
on  the  coordinate  fixed  at  the  middle  of  the  lower  surface  of  the  moving  bubble, 
are  shown  in  Figure  5.4.  The  recirculation  inside  the  bubble,  which  starts  a  short 
distance  from  the  interface,  is  observed  even  for  the  bubble  undergoing  phase  change. 
In  the  vicinity  of  the  interface,  the  streamlines  emanating  from  the  interface  are  a 
result  of  evaporation.  The  detached  recirculating  wake  is  also  observed  for  the  case 
with  higher  Reynolds  number. 

The  corresponding  temperature  profiles  around  the  bubble  for  these  two  cases 
are  shown  in  Figure  5.5.  The  thermal  boundary  layer  is  formed  where  the  boundary 
layer  is  thinner  around  the  upper  surface,  and  thicker  around  the  lower  surface  of 
the  rising  bubble  owing  to  the  relative  motion  between  the  rising  bubble  and  the 
surrounding  liquid.  It  is  also  clear  that  the  higher  the  Reynolds  number,  the  thinner 
the  thermal  boundary  layer.  The  tail-shaped  structure  in  the  wake  results  from  the 
separation  of  the  boundary  layer. 

The  local  Nusselt  numbers,  i.e.,  dimensionless  heat  flux,  along  the  bubble  surface 
at  the  end  of  the  simulations  are  depicted  in  Figure  5.6.  The  marching  direction  in 
the  figure  is  from  the  middle  of  the  lower  surface  to  the  middle  of  the  upper  surface 
of  the  bubble  for  the  Nusselt  numbers.  The  highest  heat  flux  occurs  at  certain  point 
where  the  relative  velocity  between  the  liquid  and  the  bubble  is  the  highest.   The 
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Figure  5.3:  Growth  rate  of  the  bubble  radius  for  a  rising  bubble  obtained  in  the 
present  simulation  for  cases  (a)  Re  =  10.0,  Fr  =  1.0,  We  =  0.2,  J  a  =  1.0,  Pr  =  1.0, 
Pe  =  10.0,  Pi/pv  =  1605,  fu/fr  =  22,  (b)  Re  =  174.81,  Fr  =  1.0,  We  =  0.42, 
Ja  a  16.43,  Pr  =  8.37,  Pe  =  1463.7,  p,/p„  =  527,  m/fxv  =  41. 
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Figure  5.4:  Flow  structures  at  the  end  of  the  simulations  for  cases  (a)  Re  =  10.0, 
Fr  =  1.0,  We  =  0.2,  Ja  =  1.0,  Pr  =  1.0,  Pe  =  10.0,  pt/pv  =  1605,  //,///„  =  22, 
(b)  Re  =  174.81,  Fr  =  1.0,  We  =  0.42,  Ja  =  16.43,  Pr  =  8.37,  Pe  -  1463.7, 
Pi/Pv  —  527,  \iij\iv  —  41.  The  streamlines  are  observed  on  the  reference  frame 
attached  to  the  moving  bubble. 
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spots  of  the  maximum  heat  flux  are  not  at  the  middle  of  the  upper  surface,  but 
somewhere  downstream. 

The  overall  Nusselt  number  versus  time  for  the  two  cases  are  plotted  in  Figure 
5.7.  The  average  heat  flux  almost  became  constant  which  means  the  thermal  bound- 
ary layers  are  fully  developed  at  the  time  the  simulations  ceased.  The  case  with  the 
higher  Reynolds  number  and  Jakob  number  has  higher  average  heat  flux  than  the 
other. 

The  time  evolutions  of  the  rising  velocity  of  the  bubble,  which  is  the  mean  of  the 
velocity  values  at  the  middle  of  the  lower  and  upper  surfaces  of  the  moving  bubble, 
are  shown  in  Figure  5.8.  The  rising  velocity  of  the  bubble  is  increasing  owing  to 
its  increased  volume.  It  is  noted  that  the  present  simulation  computed  the  solution 
under  the  varying  bubble  rising  velocity  resulted  from  the  buoyancy  effect.  Other 
investigations  on  the  bubble  growth  often  employ  an  uniform  mean  flow  to  simulate 
the  effect  of  the  relative  velocity  on  the  bubble  growth.  The  present  simulation  is 
closer  to  the  realistic  situations. 

In  Chapter  4,  we  studied  the  effect  of  density  ratio  on  the  bubble  dynamics  for 
Re  =  100,  andWe  =  4.  Here  we  compute  the  bubble  growth  under  the  same  Reynolds 
and  Weber  numbers. 

Figure  5.9  shows  the  development  of  bubble  shapes,  under  the  influence  of  phase 
change  and  buoyancy,  at  selected  time  instants  for  three  cases.  In  these  cases,  the  Re, 
We  and  viscosity  ratio  are  fixed,  while  the  density  ratio  is  varied  from  0.1  to  0.001. 
It  is  noted  that  the  density  ratio  directly  influences  the  Ja  and  the  interface  speed, 
as  indicated  in  (1.25).  Between  the  density  ratios  of  0.1  and  0.01,  while  the  bubble 
size  grows  slightly  faster  with  a  smaller  density  ratio,  the  shapes  are  similar  between 
the  two  cases.  However,  as  the  density  ratio  is  reduced  to  0.001,  which  is  closer  to 
a  normal  boiling  heat  transfer  case,  significant  differences  in  bubble  growth  rate  and 
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Figure  5.5:  The  temperature  profile  at  the  end  of  the  simulations  for  cases  (a)  Re  = 
10.0,  Fr  =  1.0,  We  =  0.2,  Ja  =  1.0,  Pr  =  1.0,  Pe  =  10.0,  pt/pv  =  1605,  in) pv  =  22, 
(b)  Re  =  174.81,  Fr  =  1.0,  We  =  0.42,  Ja  =  16.43,  Pr  =  8.37,  Pe  =  1463.7, 
Pi/pv  =  527,  pi/pv  =  41. 
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Figure  5.6:   The  Nusselt  number  along  the  interface  at  the  end  of  the  simulations 

for  cases  (a)  Re  =  10.0,  Fr  =  1.0,  We  =  0.2,  Ja  =  1.0,  Pr  =  1.0,  Pe  =  10.0, 

pi/pv  =  1605,  pi/pv  =  22,  (b)  Re  =  174.81,  Fr  =  1.0,  We  =  0.42,  Ja  =  16.43, 

Pr  =  8.37,  Pe  =  1463.7,  p//p„  =  527,  pt/pv  =  41. 
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Figure  5.7:  The  overall  Nusselt  number  versus  time  for  cases  (a)  Re  =  10.0,  Fr  =  1.0, 
We  =  0.2,  Ja  =  1.0,  Pr  =  1.0,  Pe  =  10.0,  pi/pv  =  1605,  m/nv  =  22,  (b)  Re  = 
174.81,  Fr  =  1.0,  We  =  0.42,  Ja  =  16.43,  Pr  -  8.37,  Pe  =  1463.7,  #//>„  =  527, 

flt/Hv  a  41. 
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(a) 


(b) 


Figure  5.8:  The  bubble  rising  velocity  versus  time  for  cases  (a)  Re  —  10.0,  Fr  =  1.0, 
We  =  0.2,  Ja  =  1.0,  Pr  =  1.0,  Pe  =  10.0,  Pl/pv  =  1605,  /*,///„  =  22,  (b)  Re  = 
174.81,  Fr  =  1.0,  We  =  0.42,  Ja  =  16.43,  Pr  =  8.37,  Pe  =  1463.7,  Pl/pv  =  527, 
Hi/Hr,  =  41. 
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Figure  5.9:   Development  of  phase  change  shapes  for  Re  =  100,   We  =  4,   Fr  — 
1,  fjLv/fj,i  =  1  with  different  density  ratios. 


interface  shape  are  observed.  Figure  5.10  shows  the  flow  structures  at  the  final  stage 
of  each  case. 
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Figure  5.10:  Flow  structures  at  the  end  of  the  simulaions  for  cases  with  Re  = 
100,  We  =  4,  Fr  =  1,  //„//*,  =  1,  Pv/Pl  =  0.1,  0.01,  0.001;  The  streamlines 
are  observed  on  the  reference  frame  attached  to  the  moving  bubble. 


CHAPTER  6 

BUBBLE  RISE  AND 
COLLAPSE  DUE  TO  CONDENSATION 

The  collapse  process  is  more  challenging  numerically  than  the  growth  case  be- 
cause of  the  fact  that  the  ratio  of  bubble  surface  area  and  its  volume  is  increasing 
with  time  in  contrast  to  decreasing  in  the  growth  case.  This  trend  in  geometrical 
relation  indicates  that  the  radius  variation  induced  by  interfacial  heat  flux  would  be 
increasing  as  the  bubble  size  is  reduced. 

In  the  growth  situation,  the  ratio  between  the  thermal  boundary  layer  thickness 
and  bubble  radius  will  eventually  reach  a  constant  value  because  they  both  grow  in 
the  same  direction.  The  thin  boundary  layer  theory  could  thus  be  used  to  find  the 
heat  flux  on  the  bubble  surface.  However  that  is  not  the  case  for  bubble  collapse. 

Since  the  thermal  boundary  layer  surrounding  a  collapsing  bubble  becomes 
thicker  as  the  bubble  size  is  reduced,  the  thin  boundary  layer  approximation  is  not 
valid  for  describing  the  heat  transfer  in  this  expanding  thermal  region. 

Unlike  bubble  growth,  there  is  no  analytical  solution  valid  for  the  entire  collaps- 
ing process.  Some  used  the  thin  boundary  layer  assumption,  such  as  Florschuetz  & 
Chao  (1965),  to  obtain  an  analytical  result  which  is  valid  only  at  the  beginning  of  the 
collapse  for  large  Jakob  numbers.  Okhotsimskii  (1988)  solved  numerically  the  heat 
transfer-controlled  collapse  based  on  the  same  theoretical  model  of  Plesset  &  Zwick 
(1954)  and  tabulated  the  time  evolution  of  the  bubble  radius  over  a  wide  range  of 
Jakob  numbers  (  10-2  <  J  a  <  103  ).  In  his  analysis,  he  isolated  Jakob  number  to  be 
the  only  dominant  factor  of  the  collapse  process. 
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For  a  translating  and  collapsing  bubble,  Tokada  et  al.  (1968),  Akiyama  (1973), 
Moalem  k  Sideman  (1973),  and  Dimic  (1977)  etc.  obtained  equations  for  predicting 
the  temporal  decrease  of  the  condensing  bubble  diameter.  They  are  reviewed  in  Chen 
&  Mayinger  (1992)  along  with  some  other  formulas. 

Chen  &  Mayinger  (1992)  gave  an  empirical  correlation  describing  the  temporal 
decrease  in  the  bubble  diameter  using  their  experimental  data  as 

7  =  (1  -  0.5QRe°-7Pr°-5JaFo)°-9  (6.1) 

in  which  7  —  R/Ro,  Fourier  number  Fo  =  ait/(2R0)2  and  the  Reynolds  number 
based  on  the  initial  bubble  diameter. 

For  the  rising  and  collapsing  bubble,  the  same  velocity  scale  uT  =  yfgL  in  the 
translating  bubble  growth  cases  is  employed.  Thereby  the  Froude  number  is  always 
1.0. 

Two  cases  were  computed  for  the  rising  and  collapsing  bubble  based  on  thermal 
properties  of  water  and  ethanol,  respectively.  The  dimensionless  parameters  for  the 
two  cases  are  summarized  in  the  following. 

1.  The  first  case  we  simulated  is  based  on  the  thermal  properties  of  Ethanol  un- 
der the  atmospheric  conditions,  as  those  listed  in  Table  1.4.  A  AT  =  10°C 
superheat  is  considered.  Under  this  superheat,  the  bubble  departure  diameter 
is  about  1.0  mm  according  to  empirical  calculations  (see  Carey,  1992),  which 
is  the  length  scale  L.  The  velocity  scale  is  thus  ur  =  y/gL  =  9.9  x  10~2  m/s. 
With  these  reference  scales  and  material  properties,  the  resulted  dimensionless 
parameters  are 

Re  =  174.81 

Fr  =  1.0 

Pe  =  1463.7 
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We  =  0.42 
Pr  =  8.37 
Ja  =  16.43 

2.  The  second  case  we  simulated  is  based  on  the  thermal  properties  of  water  under 
the  atmospheric  conditions,  as  those  listed  in  Table  1.3.  A  AT  =  3°C  superheat 
is  considered.  Under  this  superheat,  the  bubble  departure  diameter  is  about 
0.9  mm  according  to  empirical  calculations  (see  Carey,  1992),  which  is  used  as 
the  length  scale  L.  The  velocity  scale  is  thus  ur  =  yfgL  =  9.39  x  10-2  m/s. 
With  these  reference  scales,  the  resulted  dimensionless  parameters  are 

Re  =  291.86 

Fr  =  1.0 

Pe  =  503.41 

We  =  0.13 

Pr  =  1.72 

Ja  =  9.0 

In  Figure  6.1,  the  temporal  variation  of  bubble  radius  for  case  1  obtained  by 
the  present  simulation,  correlation  of  (6.1)  and  prediction  of  Akiyama  (1973)  are 
all  plotted.  In  our  simulation,  we  compute  the  flow  fields  both  inside  and  outside 
the  bubble.  A  certain  number  of  grid  cells  have  to  be  enclosed  inside  the  bubble. 
Therefore  we  were  unable  to  compute  the  entire  collapse  process  but  to  stop  at  a 
certain  size  to  ensure  that  we  have  enough  cells  to  resolve  the  flow  fields  inside  the 
bubble.  Thereby  the  dimensionless  bubble  radius  didn't  reach  zero  at  the  end  of 
curves. 

Figure  6.1  indicates  that  the  bubble  collapse  speed  is  faster  in  the  earlier  stage 
than  that  predicted  by  Chen  &  Mayinger  (1992).   After  that,  the  collapse  speed  is 
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Figure  6.1:  Comparison  of  the  temporal  change  of  bubble  radius  predicted  by  (6.1), 
Akiyama  (1973)  and  the  present  simulation  for  Re  =  174.81,  Fr  =  1.0,  We  =  0.42, 
Ja  =  16.43,  Pr  =  8.37,  Pe  =  1463.7,  pt/pv  =  527,  m/fa  =  41;  7  =  R/R0,  Fo  = 
ait/(2R0)2. 
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Figure  6.2:  Comparison  of  the  temporal  change  of  bubble  radius  predicted  by  (6.1), 
Florschuetz  &  Chao  (1965),  Akiyama  (1973)  and  the  present  simulation  for  Re  = 
291.86,  Fr  =  1.0,  We  =  0.13,  Ja  =  9.0,  Pr  =  1.72,  Pe  =  503.41,  pt/pv  =  1605, 
m/fr  =  22;  7  =  R/Ro,  Fo  =  att/(2R0)2. 


slower  during  the  majority  of  the  collapse  process  than  that  of  Chen  &  Mayinger 
(1992).  Moreover,  present  results  are  closer  to  those  by  Chen  &  Mayinger  (1992) 
which  were  obtained  from  measurements.  Akiyama  (1973)  used  the  heat  transfer 
rate  similar  to  that  around  a  moving  solid  sphere  to  obtain  his  analytical  solution, 
which  probably  resulted  in  a  slower  collapse  rate  than  that  of  a  deformed  bubble. 

In  Figure  6.2,  the  temporal  variation  of  bubble  radius  for  case  2  obtained  by  the 
present  simulation,  correlation  of  (6.1)  and  prediction  of  Akiyama  (1973)  as  well  as 
Florschuetz  &  Chao  (1965)  are  plotted. 
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From  these  two  cases,  the  observation  is  that  the  present  simulation  is  not  in 
favorable  agreement  with  the  correlation  proposed  by  Chen  &  Mayinger  (1992)  for 
the  time  evolution  of  the  bubble  radius  of  a  collapsing  and  rising  bubble. 

Chen  &  Mayinger  (1992)  also  suggested  that  their  data  were  better  fitted  by 
(6.1)  for  values  of  Pe/Ja  equal  to  several  hundred.  When  the  ratio  of  Pe/Ja  drops 
below  60,  a  significant  disagreement  is  observed  (see  Legendre,  Boree  &;  Magnaudet, 
1998).  This  is  indeed  the  outcome  from  our  simulations,  case  1  has  a  Pe/Ja  ratio 
over  60  so  that  a  better  match  is  observed.  Case  2  has  a  ratio  lower  than  60,  and 
shows  substantial  disagreement. 


CHAPTER  7 
SUMMARY  AND  FUTURE  WORK 

7.1     Summary  of  Present  Work 

The  objective  of  the  present  work  has  been  to  develop  a  predictive  numerical 
capability  for  single  bubble  dynamics  with  phase  change.  In  order  to  do  this,  the 
fixed-grid,  sharp  interface  methodology  has  been  developed  (Ye  et  a/.,  20016, a).  The 
major  components  developed  in  this  work  for  the  entire  numerical  simulation  method 
are 

•  AC2  cubic  spline  curve  fitting  in  conjunction  with  the  fairing  algorithm  is  im- 
plemented to  numerically  represent,  track  the  interface,  and  obtain  the  accurate 
geometric  informations.  The  critical  advantage  of  this  scheme  is,  by  eliminat- 
ing the  oscillations  in  the  curvature  calculation,  to  enable  the  success  of  the 
computations  since  the  curvature  is  the  most  notorious  geometric  parameter 
which  often  leads  to  the  failure  of  the  overall  simulation  of  the  interfacial  flows 
in  which  the  capillary  force  is  important.  This  is  the  first  time  the  efficient  and 
capable  geometric  fairing  algorithm  is  exploited  and  implemented  to  solve  the 
curvature  calculation  problem  in  the  area  of  computational  fluid  dynamics. 

•  A  Cartesian  grid  based  flow  solver  is  developed  for  solving  the  coupled  govern- 
ing equations  on  a  fixed  grid  with  curved  boundaries.  Although  this  part  is 
developed  as  an  integrated  component  of  the  present  overall  method,  it  alone  is 
an  good  alternative  to  body-fitted  grid  approach  for  solving  flows  over  complex 
geometries.  The  major  advantage  of  this  method  is  that  it  eliminates  the  grid 
generation  phase  and  thereby  results  in  a  more  efficient  solution  procedure. 
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•  The  interface  updating  schemes  and  procedures  are  implemented  and  developed 
to  allow  the  accurate  prediction  of  the  exact  interface  shape  which  satisfies  all 
dynamical  and  thermal  conditions. 

•  All  major  components  of  this  method  are  tested  or  validated  through  compar- 
isons with  the  established  results  and  grid  resolution  study.  The  capability  of 
the  method  is  demonstrated  in  simulations  of  buoyancy-driven  bubble  motion, 
bubble  growth  and  bubble  collapse  under  the  gravity.  These  are  extremely 
complex  physical  phenomena  where  multiple  interacting  effects  are  required  to 
be  modeled  simultaneously  in  addition  to  a  constantly  evolving  interface. 

7.2     Contributions  of  Present  Work 

1.  A  fixed-grid  direct  numerical  simulation  method  has  been  developed  and  verified 
for  bubble  dynamics,  heat  transfer  and  phase  change  phenomena.  The  original 
contribution  of  this  research  is  the  first  introduction  of  a  sharp  interface  between 
the  vapor  and  liquid  phases,  which  allows  the  simulation  of  density  ratios  up 
to  the  order  of  1000.  The  interface  location  is  obtained  as  part  of  the  solution 
which  accounts  for  all  coupled  dynamic  and  thermal  effects. 

2.  The  current  method  is  believed  to  be  the  only  method  available  at  present 
which  resolves  the  curvature  calculation  issue,  treats  the  interface  as  a  sharp 
discontinuity,  honors  the  mass  conservation  principle,  is  capable  of  handling  the 
large  property  ratios  and  phase  change. 

3.  For  the  isothermal  cases,  the  predictions  of  drag  coefficients  and  the  shape 
deformation  of  bubbles  for  Reynolds  numbers  ranging  from  2  to  100  and  We- 
ber numbers  from  2  to  20  were  shown  independently  to  agree  with  previous 
published  results  of  Ryskin  &  Leal  (1984). 
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4.  It  was  also  demonstrated  that  the  current  method  is  capable  of  predicting 
accurate  results  for  a  wide  range  of  Re,  We,  Ja,  Pe  numbers  and  property 
jumps  at  the  interface. 

5.  For  a  stationary  and  growing  bubble,  the  predicted  growth  rate  was  found  to 
agree  with  the  theoretical  limit  of  oc  i1/2  (see  Appendix  A).  However,  for  a  rising 
and  growing  bubble,  the  predicted  growth  rate  was  approaching  the  theoretical 
limit  of  a  t2?3  (see  Appendix  A)  for  a  spherical  bubble  but  did  not  reach  it 
exactly  owing  to  bubble  deformation. 

6.  The  empirical  correlation  proposed  for  bubble  collapse  is  assessed  by  the  sim- 
ulations in  the  present  work.  The  present  numerical  model  predicts  that  the 
dynamics  of  a  collapsing  bubble  is  qualitatively,  but  not  quantitatively  similar 
to  the  empirical  correlation  proposed  by  Chen  &  Mayinger  (1992). 

7.3     Future  Work 
7.3.1     Applications 

There  are  many  important  phenomena  in  our  nature  can  be  modeled  as  multi- 
phase, multicomponent  flows  with  dynamic  interfaces  such  as  fluid-fluid,  fluid-solid 
interaction  problems.  Compared  to  pure,  single  phase  materials  and  fluids,  the  nu- 
merical simulations  of  multiphase  flows  with  dynamic  interfaces  are  among  the  most 
difficult  computational  problems  because  of  the  major  challenge  of  resolving  the  in- 
terfacial  motion.  In  this  research,  a  numerical  capability  has  been  developed  to  begin 
probing  this  frontier.  The  method  is  applied  to  study  the  vapor  bubble  dynamics 
which  is  just  one  of  examples  of  flows  with  dynamic  interfaces. 

A  target  application  of  the  present  method  is  to  investigate  the  flows  in  the  bio- 
logical systems.  In  biomedical  engineering,  investigations  of  the  interaction  of  blood 
cells  with  elastic  vessels  are  underway  to  explore  the  mechanical  mechanisms  of  heart 
and  vessel  diseases.  The  study  of  blood  flow  in  the  small  artificial  organs  helps  develop 
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better  devices  with  reduced  blood  damage  when  blood  circulates  through  them.  For 
these  types  of  problems,  blood  can  be  modeled  as  numerous  deformable  cells  (largely 
red  blood  cells)  suspended  in  a  Newtonian  fluid  (plasma).  The  blood  cells  are  bod- 
ies containing  fluid  gel  within  a  solid  membrane.  The  blood  damage  occurs  when 
those  blood  cells  are  broken  due  to  unsustainable  stress  and  deformation  caused  by 
interaction  among  the  blood  cells,  plasma  and  the  solid  walls.  The  dynamics  of  the 
blood  cell  membranes  are  central  to  the  cell  damage  process.  In  order  to  quantify  the 
membrane  stresses  and  the  deformations,  the  microscopic  computations  at  each  cell 
level  are  necessary  to  resolve  the  numerous  dynamically  deforming  cellular  interfaces 
simultaneously. 

7.3.2     Numerical  Enhancements 

The  present  method  is  developed  to  resolve  the  interfacial  dynamics  within  a 
fixed  grid  frame.  The  sharp  interface  treatment  renders  it  well  suitable  for  simulating 
the  problems  outlined  in  the  above.  The  implementation  of  geometric  algorithms 
can  handle  numerous  interfaces  as  described  in  the  Chapter  2.  However,  several 
enhancements  have  to  be  incorporated  into  the  present  numerical  methodology. 

The  large  number  of  blood  cells  needed  to  be  covered  in  a  computation  so  as  to 
have  practically  meaningful  results,  plus  adequate  spatial  resolution  necessitate  the 
use  of  massively  parallel  computers  where  thousands  of  processors  can  each  compute 
a  sub-portion  of  the  problem.  In  this  regard,  the  present  method  has  a  significant 
advantage,  that  is,  the  parallelization  task  is  relatively  easy  owing  to  the  static-grid 
frame.  The  parallelization  issues  for  fixed  grid  solver  are  reasonably  well  understood. 

On  the  other  hand,  even  for  single  interface  problem  such  as  the  single  vapor 
bubble  dynamics  studied  in  this  thesis  research,  current  overall  solution  procedure 
relies  on  extensive  iteration  at  multiple  levels.  The  computing  cycle  often  takes  hours 
on  single  processor  computer.  The  leap  in  modeling  capability  and  quality  of  the 
present  method  is  thereby  expected  by  parallelization  of  the  numerical  algorithms. 
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The  fixed,  Cartesian  grid  can  be  readily  decomposed  into  sub-domains  which  are 
mapped  to  multiple  processors.  Each  processor  then  operates  independently  on  the 
portion  assigned  to  it,  communicating  only  the  much  smaller  shared  border  data 
when  necessary. 

The  message  passing  approach  is  the  primary  model  for  writing  parallel  pro- 
grams. The  program  execution  is  divided  among  a  number  of  processes  which  need 
not  be  running  on  different  processors.  The  process  is  a  more  general  term  referring 
to  an  independent  instruction  and  data  flow.  Each  process  thereby  can  work  on  some 
local  data.  The  data  sharing  between  processes  is  achieved  by  the  message  passing, 
that  is,  by  explicitly  sending  and  receiving  data  between  processes. 

On  modern  parallel  computers,  a  library  of  functions  or  subroutines  called  MPI 
(Message  Passing  Interface)  is  provided  for  programmer  to  insert  into  source  code 
to  perform  data  communication  between  processes.  The  programmer  must  explicitly 
divide  data  and  work  and  manage  the  communications  among  the  processes  to  achieve 
parallel  execution. 

The  merger/breakup  module  for  handling  the  interface  interaction  and  topologi- 
cal change  is  needed  for  computing  multiple  interfaces  problems.  The  merger/breakup 
function  is  implemented  in  the  present  code  and  yet  to  be  tested. 

Turbulence  might  be  important  at  high  Reynolds  number.  The  direct  numerical 
simulation  of  the  turbulence  is  still  out  of  reach  even  with  parallel  computers  and 
algorithms  because  of  the  extremely  high  spatial  resolution  requirement.  The  LES 
(Large  Eddy  Simulation)  is  a  good  compromise  which  solves  for  the  flow  structures  the 
grid  size  can  resolve,  and  models  those  smaller  than  the  grid  size.  Various  modeling 
approaches  have  been  proposed  for  those  sub-grid  eddies.  The  most  reasonable  model, 
that  is,  dynamic  sub-grid  scale  model  of  LES  has  been  implemented  in  the  research, 
yet  to  be  tested  and  incorporated  into  the  present  method. 


APPENDIX  A 
HEAT  TRANSFER-CONTROLLED  LIMITS  OF  SPHERICAL  BUBBLE  GROWTH 

For  a  spherical  bubble  undergoing  heat  transfer-controlled  growth,  the  overall 
energy  balance  for  the  evaporation  at  the  bubble  surface  could  be  expressed  as 


di 


l*R3(t) 


Xpv  =  hAT4nR2{t) 


(A.1) 


where  R(t)  is  the  instantaneous  bubble  radius,  A  is  the  latent  heat,  pv  is  the  vapor 
density,  h  is  the  heat  transfer  coefficient,  AT  is  the  degree  of  superheat. 

To  estimate  h,  Ranz  &  Marshall  (1952) 's  correlation,  considering  both  the  con- 
duction and  convection  effects,  may  be  used 
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(A.2) 


where  the  first  term  2  represents  the  contribution  from  conduction,  the  second  term 
0.6Re^Pr3  denotes  the  contribution  by  the  convection. 
From  Eq.  (A.2),  we  have 
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(A.3) 


A.l     Conduction  Dominated  Growth 


If  only  the  first  term  of  Eq.  (A.3)  is  substituted  for  h,  Eq.  (A.l)  becomes 

AkR2<^-\Pv  =  2^AT4ttR2  (A.4) 

at  R 
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Integration  of  Eq.  (A. 4)  yields 

R(t)  =  Ctl  (A.5) 

where  C  includes  all  parameters  that  are  not  function  of  the  time  such  as  thermal 
properties,  Prandtl  number  and  degree  of  superheat,  etc. 

(A.5)  indicates  that  the  conduction  dominated  growth  rate  is  proportional  to 
the  square  root  of  time. 

A. 2     Convection  Dominated  Growth 

With  only  the  second  term  of  Eq.    (A. 3)  substituted  for  h,  Eq.    (A.l)  then 
becomes 

4^2 ^A     =  0  6  («tR\  2  Prk^AT4TTR2  (A.6) 

at  \  v  )  R 

Integration  of  Eq.  (A. 4)  gives 

R(t)  =  Ct\  (A.7) 

From  (A.7),  the  convection  dominated  growth  rate  is  proportional  to  the  two- 
thirds  power  of  time. 
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